All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
BCalTagEfficiency.cc
Go to the documentation of this file.
1 //Filename RootBCalTagEfficiency.cc
2 
3 #include "BCalTagEfficiency.h"
4 #include <iostream>
5 #include <cmath>
6 #include <stdio.h>
7 
8 
9 #ifdef MARLIN_USE_AIDA
10 #include <marlin/AIDAProcessor.h>
11 #include <AIDA/IHistogramFactory.h>
12 #include <AIDA/ICloud2D.h>
13 #include <AIDA/IHistogram2D.h>
14 #endif
15 
16 #include <EVENT/LCCollection.h>
17 #include <EVENT/MCParticle.h>
18 #include <EVENT/ReconstructedParticle.h>
19 #include "IMPL/LCCollectionVec.h"
20 #include <IMPL/ReconstructedParticleImpl.h>
21 #include <IMPL/ClusterImpl.h>
22 #include <IMPL/LCRelationImpl.h>
23 #include <CLHEP/Vector/LorentzVector.h>
24 #include <EVENT/LCRelation.h>
25 #include <UTIL/LCRelationNavigator.h>
26 
27 #include <marlin/Global.h>
28 #include <gear/GEAR.h>
29 #include <gear/GearParameters.h>
30 #include <gear/BField.h>
31 #include "gear/CalorimeterParameters.h"
32 
33 // include what BcEnergyDensity class needs
34 //#include "BcEnergyDensity.h"
35 
36 // include root headers
37 #include "TRandom3.h"
38 
39 // ----- include for verbosity dependend logging ---------
40 #include "marlin/VerbosityLevels.h"
41 
42 
43 using namespace lcio ;
44 using namespace marlin ;
45 using namespace std;
46 using namespace CLHEP ;
47 
48 // Declaring the function to call FORTRAN routine
49 
50 
51 extern "C" {
52 
53  void bcalhit_(float pin_[3], float& q_,float vin_[3] ,
54  float& B_, float& eBeam_, float& zbcal_, float pout_[3], float vout_[3]);
55 }
56 
57 
58 BCalTagEfficiency aBCalTagEfficiency ;
59 
60 
61 BCalTagEfficiency::BCalTagEfficiency() : Processor("BCalTagEfficiency") {
62 
63  // modify processor description
64  _description = "ForwardVeto propagates electrons or photons from IP to the BeamCal, and it gives the Energy and Momentum of the detected ones to the output ReconstructedParticles collection" ;
65 
66 
67  // register steering parameters: name, description, class-variable,
68  // default value
69 
70  registerInputCollection( LCIO::MCPARTICLE,
71  "CollectionName" ,
72  "Name of the MCParticle collection" ,
73  _MCParticleName ,
74  std::string("MCParticlesSkimmed") ) ;
75 
76  registerInputCollection( LCIO::LCRELATION,
77  "BCALInputTruthLinkName" ,
78  "Name of the truth link input collection (from SGV or Mokka)" ,
79  _BCALInputTruthLinkName ,
80  std::string("BCALMCTruthLink") ) ;
81 
82  registerOutputCollection( LCIO::RECONSTRUCTEDPARTICLE,
83  "BCALParticleName" ,
84  "Name of the tagged e/gamma collection" ,
85  _BCALParticleName ,
86  std::string("BCALEffiParticles") ) ;
87 
88  registerOutputCollection( LCIO::CLUSTER,
89  "BCALClusterName" ,
90  "Name of the collection of clusters of tagged e's/gammas" ,
91  _BCALClusterName ,
92  std::string("BCALEffiClusters") ) ;
93 
94 
95  registerOutputCollection( LCIO::LCRELATION,
96  "BCALEffiMCTruthLinkName",
97  "Name of the tagged e/gamma to mc-particle relation collection" ,
98  _BCALEffiMCTruthLinkName,
99  std::string("BCALEffiMCTruthLink") ) ;
100 
101 
102  registerProcessorParameter("BackgroundFilename" ,
103  "Name of input file for background energy density" ,
104  backgroundfilename ,
105  std::string("bg_aver_LDC_4T_14mrad_AntiDID.root") ) ;
106 
107  // eBeam = Energy of the beam in CMS (GeV)
108  registerProcessorParameter("eBeam" ,
109  "" ,
110  eBeam,
111  float(250));
112 
113 
114  // zbcal = z point from the IP (mm)
115  registerProcessorParameter("zbcal" ,
116  "position of the BeamCal" ,
117  zbcal,
118  float(3594.9));
119 
120  registerProcessorParameter("thresholdmin",
121  "Minimum energy of detected particles",
122  thresholdMin,
123  float(40));
124 
125 
126  registerProcessorParameter("thresholdmax",
127  "Maximum energy of detected particles" ,
128  thresholdMax,
129  float(260));
130 
131 
132  registerProcessorParameter("smearEnergy",
133  "true: electrons/photons written to new collection with smeared energy (and momentum)",
134  smearEnergy,
135  bool(1));
136 
137  registerProcessorParameter("detectAll",
138  "true: all electrons/photons in BCAL acceptance will be written to new collection",
139  detectAll,
140  bool(1));
141 
142  registerProcessorParameter("densityScaling",
143  "Scaling of the energy density w.r.t default setup",
144  densityScaling,
145  float(1.));
146 
147  registerProcessorParameter("DBDsample",
148  "true: assume DBD-style event sample with MCParticles with crossing angle, false: LoI-Style, MCParticles head-on",
149  DBDsample,
150  bool(1));
151 
152  registerProcessorParameter("newMap",
153  "true: new map with bcal layers 1 - 30, pair mon layer 31; false: old map with pair monitor in layer1",
154  newMap,
155  bool(1));
156 
157  registerProcessorParameter("writeTree",
158  "true: write root tree with debug info",
159  writeTree,
160  bool(1));
161 
162  registerProcessorParameter("writeSGVMap",
163  "true: write ascii map for SGV",
164  writeSGVMap,
165  bool(0));
166 
167  registerProcessorParameter("SGVMapFilename",
168  "name of file for SGV map",
169  SGVmapfilename,
170  std::string("SGVmap.root") );
171 
172  registerProcessorParameter("useBCALInputTruthLink",
173  "true: read BCALMCTruthLink from SGV / Mokka to find correct MCParticles",
174  useInputClusters,
175  bool(0));
176 
177 }
178 
179 void BCalTagEfficiency::init() {
180 
181  streamlog_out(DEBUG) << " init called "
182  << std::endl ;
183 
184  // instanciate background handler object
185  bc_en = new BcEnergyDensity(backgroundfilename.c_str());
186 
187  streamlog_out(DEBUG) << " background initialized "
188  << std::endl ;
189 
190  bField = Global::GEAR->getBField().at(gear::Vector3D(0.,0.,0.)).z();
191  streamlog_out(MESSAGE) << "BCalTagEfficiency::init: B-FIELD = " << bField << std::endl;
192 
193  //read parameters from gear file
194  const gear::CalorimeterParameters& bcparam = Global::GEAR->getBeamCalParameters(); //in bcparam are read all parameters from the gearfile; the gear file is already opened by default by Marlin (the name is known from the Marlin steering file)
195  xingangle = bcparam.getDoubleVal("beam_crossing_angle");
196 
197  // parameters of the Lorentz transformation / rotation matrix
198  alpha = xingangle/2000; // half x-ing angle in radians
199  gamma = sqrt(1 + tan(alpha)*tan(alpha));
200  betagamma = tan(alpha);
201 
202  // usually a good idea to
203  printParameters() ;
204 
205  _nRun = 0 ;
206  _nEvt = 0 ;
207 
208  PI = 4.*atan(1.);
209 
210  if (writeTree) {
211 
212  // root stuff
213  rootfile = new TFile("BCalEffi.root","recreate");
214  tree = new TTree("Treename","Tree for electrons");
215 
216  tree->Branch("mcp", &mcp,"mcp/I");
217  tree->Branch("pdg", pdg,"pdg[mcp]/I");
218  tree->Branch("energy", energy,"energy[mcp]/F");
219  tree->Branch("ePrime", ePrime,"ePrime[mcp]/F");
220  tree->Branch("pxPrime",pxPrime,"pxPrime[mcp]/F");
221  tree->Branch("pxIP", pxIP,"pxIP[mcp]/F");
222  tree->Branch("pyIP", pyIP,"pyIP[mcp]/F");
223  tree->Branch("pzIP", pzIP,"pzIP[mcp]/F");
224  tree->Branch("scaleP", scaleP,"scaleP[mcp]/F");
225  tree->Branch("phiIP", phiIP,"phiIP[mcp]/F");
226  tree->Branch("theIP", theIP,"theIP[mcp]/F");
227  tree->Branch("lposx", lposx,"lposx[mcp]/F");
228  tree->Branch("lposy", lposy,"lposy[mcp]/F");
229  tree->Branch("lposz", lposz,"lposz[mcp]/F");
230  tree->Branch("gposx", gposx,"gposx[mcp]/F");
231  tree->Branch("gposy", gposy,"gposy[mcp]/F");
232  tree->Branch("gposz", gposz,"gposz[mcp]/F");
233  tree->Branch("radius", radius,"radius[mcp]/F");
234  tree->Branch("phi", phi,"phi[mcp]/F");
235  tree->Branch("ebkg", ebkg,"ebkg[mcp]/F");
236  tree->Branch("ebkg_err",ebkg_err,"ebkg_err[mcp]/F");
237  tree->Branch("efficiency", efficiency,"efficiency[mcp]/F");
238  tree->Branch("rand" , rand,"rand[mcp]/F");
239  tree->Branch("tag" , tag,"tag[mcp]/I");
240  }
241 
242  if (writeSGVMap) {
243  nbinx = 50;
244  xmin = -140.;
245  xmax = 160.;
246  nbiny = 50;
247  ymin = -150.;
248  ymax = 150.;
249  TFile *SGVmaprootfile = new TFile(SGVmapfilename.c_str(),"recreate");
250  SGVmapP = new TH2D("hmapP","mapP", nbinx, xmin, xmax, nbiny, ymin, ymax);
251  SGVmapN = new TH2D("hmapN","mapN", nbinx, xmin, xmax, nbiny, ymin, ymax);
252  SGVmapP->Sumw2();
253  SGVmapN->Sumw2();
254  double wx = (xmax-xmin)/nbinx;
255  double wy = (ymax-ymin)/nbiny;
256  double vloc[3]; // local coordinates
257  for (int iz = -1; iz < 2; iz += 2) {
258  float halfxingangle = iz * alpha;
259  for (int ix = 0; ix < nbinx; ++ix) {
260  // bin center in x
261  double x = xmin + (ix+0.5)*wx;
262  // convert to local: rotation around y axis of -7 mrad.
263  // rotation goes in different directions for +z / -z
264  vloc[0] = cos(-halfxingangle)* x + sin(-halfxingangle)* (iz*zbcal);
265  vloc[2] = -sin(-halfxingangle)* x + cos(-halfxingangle)* (iz*zbcal);
266  for (int iy = 0; iy < nbiny; ++iy) {
267  // bin center in y
268  vloc[1] = ymin + (iy+0.5)*wy;
269  // Change from cartesian to cylindrical coordinates !LOCAL!
270  Hep3Vector lvec (vloc[0], vloc[1], vloc[2]);
271  double rloc = lvec.perp();
272  double philoc = lvec.phi(); // gives [-PI, PI]
273  if (philoc < 0) philoc += 2*PI; // have [0, 2 PI] now
274 
275  double sumEDens = 0;
276  double sumEDensErr = 0;
277  double EDens = 0;
278  double EDensErr = 0;
279  bool ok = false;
280 
281  for (int layer=1; layer<31; layer++) {
282  if (!newMap && layer ==1) continue; // skip layer if using old map!
283  ok = bc_en->GetEnergyDensity( layer * iz, rloc, philoc, &EDens, &EDensErr);
284  if (ok) {
285  sumEDens += EDens;
286  sumEDensErr += EDensErr/55.*EDensErr/55.;
287  }
288  // in keyhole region
289  else {
290  sumEDens = -1.;
291  sumEDensErr = -1.;
292  layer = 31;
293  }
294  }
295  sumEDensErr = (sumEDensErr>=0) ? sqrt(sumEDensErr) : -1.;
296 
297 
298  //Scaling of the Energy density with scaling factor
299  if (sumEDens > 0) {
300  sumEDens *= densityScaling;
301  sumEDensErr *= densityScaling;
302  }
303  if (iz > 0) SGVmapP->Fill(x,vloc[1],sumEDens);
304  else SGVmapN->Fill(x,vloc[1],sumEDens);
305  streamlog_out(MESSAGE) << "SGVmap " << x << " " << vloc[1] << " " << iz*zbcal << " " << sumEDens << " " << sumEDensErr << std::endl;
306  }
307  }
308  }
309 
310  SGVmapP->Write();
311  SGVmapN->Write();
312  SGVmaprootfile->Write();
313  SGVmaprootfile->Close();
314  //delete SGVmaprootfile;
315  //delete SGVmapP;
316  //delete SGVmapN;
317 
318  } // if write SGV map
319 
320 }
321 
322 void BCalTagEfficiency::processRunHeader( LCRunHeader* run) {
323 
324  _nRun++ ;
325 }
326 
327 void BCalTagEfficiency::processEvent( LCEvent * evt ) {
328 
329  streamlog_out(DEBUG) << "+++++++++++++++++ start processing event: " << evt->getEventNumber()
330  << " in run: " << evt->getRunNumber()
331  << std::endl ;
332 
333  // Create output collection
334  LCCollectionVec* BCALCol = new LCCollectionVec(LCIO::RECONSTRUCTEDPARTICLE);
335  //LCCollectionVec* BCALEffiMCTruthLink = new LCCollectionVec(LCIO::LCRELATION);
336  LCCollection* BCALEffiMCTruthLink = 0;
337  LCCollectionVec* BCALClusters = new LCCollectionVec(LCIO::CLUSTER);
338  LCRelationNavigator bcalTruthRelNav(LCIO::RECONSTRUCTEDPARTICLE , LCIO::MCPARTICLE ) ;
339 
340  LCCollection* col = 0;
341  // read MCParticles and track to BeamCal
342  if (!useInputClusters) {
343  streamlog_out(DEBUG) << "Using MCParticle collection!" << std::endl;
344  col = evt->getCollection( _MCParticleName ) ;
345  if (!col) streamlog_out(ERROR) << "MCParticle collection not found!!!" << std::endl;
346  }
347  // loop over BCalParticles, find corresponding MCParticle, track to BeamCal
348  else {
349  streamlog_out(DEBUG) << "Using BCALTruthLink collection!" << std::endl;
350  col = evt->getCollection( _BCALInputTruthLinkName ) ;
351  if (!col) streamlog_out(ERROR) << "TruthLinker not found!!!" << std::endl;
352  }
353 
354  if( col != 0 ){
355 
356  int nMCP = col->getNumberOfElements() ;
357 
358  // Arguments to be called from the FORTRAN routine
359  float q =0;
360  float pin[3]= {0, 0, 0}; // momenta at IP, !in coord sys of input!
361  float vin[3]= {0, 0, 0}; // coordinates of PCA, !in coord sys of input!
362  float gvout[3]= {0, 0, 0}; // coordinates at BeamCal surface, global coord sys
363  float lvout[3]= {0, 0, 0}; // coordinates at BeamCal surface, local coord sys
364  float pout[3]= {0, 0, 0}; // momenta at BeamCal surface !in coord sys of input! WILL NOT BE USED!
365 
366  mcp = 0;
367 
368  for(int i=0; i < nMCP ; i++){
369 
370  MCParticle* p = 0;
371  LCRelation* rp = 0;
372 
373  if (!useInputClusters) {
374  // take all genStat = 1 electrons & photons
375  p = dynamic_cast<MCParticle*>( col->getElementAt( i ) ) ;
376  if ( p->getGeneratorStatus()!=1 ) continue;
377  pdg[mcp] = p->getPDG();
378  if ( abs(pdg[mcp]) != 11 && abs(pdg[mcp]) != 22 ) continue;
379  streamlog_out(DEBUG) << "found MCP..." << std::endl;
380  }
381  else {
382  streamlog_out(DEBUG) << "looking for MCP..." << std::endl;
383  // pick out correct MCP using SGV / Mokka information!
384  rp = dynamic_cast<LCRelation*>( col->getElementAt( i ) ) ;
385  if (!rp) {
386  streamlog_out(WARNING) << " entry in LCRelation not found! " << std::endl;
387  continue;
388  }
389  if (rp->getWeight()<0.5) {
390  streamlog_out(WARNING) << " relation weight low: " << rp->getWeight() << std::endl;
391  continue;
392  }
393  p = dynamic_cast<MCParticle*>(rp->getTo());
394  streamlog_out(DEBUG) << "found MCP with PDG = " << p->getPDG() ;
395  streamlog_out(DEBUG) << " and energy = " << p->getEnergy() << std::endl;
396  }
397 
398  // Get initial energy
399  energy[mcp] = p->getEnergy();
400 
401  // Only use the electron if it's energy > 40 GeV
402  if((energy[mcp] > thresholdMin && energy[mcp] < thresholdMax)) {
403 
404  // get momentum
405  pin[0] = p->getMomentum()[0];
406  pin[1] = p->getMomentum()[1];
407  pin[2] = p->getMomentum()[2];
408 
409  pxIP[mcp] = pin[0];
410  pyIP[mcp] = pin[1];
411  pzIP[mcp] = pin[2];
412 
413  // this is with / without x-ing angle for DBD / LoI samples
414  HepLorentzVector p4v(pin[0], pin[1], pin[2], energy[mcp]);
415  phiIP[mcp] = p4v.phi();
416  streamlog_out(DEBUG) << "phi of electron/photon before +=2PI " << phiIP[mcp] << std::endl;
417  if (phiIP[mcp] < 0) phiIP[mcp] += 2*PI;
418  theIP[mcp] = p4v.theta();
419 
420  // theta larger than 0.06 (r/z = 200mm / 3500m = 0.057) will go into main detector
421  if (theIP[mcp] > 0.06 && theIP[mcp] < PI-0.06 || theIP[mcp] < 0.003 || theIP[mcp] > PI-0.003) continue;
422  streamlog_out(DEBUG) << "energy of electron/photon = " << energy[mcp] << std::endl;
423  streamlog_out(DEBUG) << "phi at IP of electron/photon = " << phiIP[mcp] << std::endl;
424  streamlog_out(DEBUG) << "pt of electron/photon = " << p4v.perp() << std::endl;
425 
426  // gets de particle charge
427  q = p->getCharge();
428 
429  // call to the FORTRAN routine bcalhit
430  // gives pout and vout
431 
432  if( abs(pdg[mcp]) == 11) { // electron, full B field tracking
433  bcalhit_(pin,q,vin,bField,eBeam,zbcal,pout,gvout);
434  }
435  else { // photon - why do we have to treat this separately?
436  gvout[0]=(pin[0]/abs(pin[2]))*zbcal;
437  gvout[1]=(pin[1]/abs(pin[2]))*zbcal;
438  gvout[2]=(pin[2]/abs(pin[2]))*zbcal;
439  pout[0]=pin[0];
440  pout[1]=pin[1];
441  pout[2]=pin[2];
442  }
443 
444  // rotation goes in different directions for +z / -z
445  float halfxingangle = pin[2]>0 ? alpha : -alpha;
446 
447  // now sort out whether we have crossing angle already in MCParticles or not
448  if (DBDsample) {
449  // calculate local BeamCal coordinates in order to compare with background map
450  // Rotation around y axis of -7 mrad. (back from global to local!)
451  lvout[0] = cos(-halfxingangle)* gvout[0] + sin(-halfxingangle)* gvout[2];
452  lvout[1] = gvout[1];
453  lvout[2] = -sin(-halfxingangle)* gvout[0] + cos(-halfxingangle)* gvout[2];
454 
455  // don't need lpout? -> check!
456  }
457  else {
458  // move gvout to local variables
459  // (global coordinates without crossing-angle correspond directly to local
460  // BeamCal coordinates!)
461  // calculate global coordinates
462 
463  lvout[0] = gvout[0];
464  lvout[1] = gvout[1];
465  lvout[2] = gvout[2];
466 
467  // Rotation around y axis of 7 mrad. (forward from local to global!)
468  gvout[0] = cos(halfxingangle)* lvout[0] + sin(halfxingangle)* lvout[2];
469  gvout[1] = lvout[1];
470  gvout[2] = -sin(halfxingangle)* lvout[0] + cos(halfxingangle)* lvout[2];
471 
472  }
473 
474  streamlog_out(DEBUG) << "entry point on BeamCal surface, global ILD coordinates:" << std::endl;
475  streamlog_out(DEBUG) << "gvout[0] = " << gvout[0] << ", gvout[1] = " << gvout[1] << ", gvout[2] = " << gvout[2] << std::endl;
476  streamlog_out(DEBUG) << "entry point on BeamCal surface, local coordinates:" << std::endl;
477  streamlog_out(DEBUG) << "lvout[0] = " << lvout[0] << ", lvout[1] = " << lvout[1] << ", lvout[2] = " << lvout[2] << std::endl;
478 
479 
480  // Change from cartesian to cylindrical coordinates !LOCAL!
481  Hep3Vector lvvec (lvout[0], lvout[1], lvout[2]);
482 
483  radius[mcp] = lvvec.perp();
484  streamlog_out(DEBUG) << "radius on BeamCal surface in local coordinates = " << radius[mcp] << std::endl;
485 
486  phi[mcp] = lvvec.phi(); // gives [-PI, PI]
487  streamlog_out(DEBUG) << "phi on BeamCal surface in local coordinates = " << phi[mcp] << std::endl;
488  if (phi[mcp] < 0) phi[mcp] += 2*PI; // have [0, 2 PI] now
489  streamlog_out(DEBUG) << "phi after +2PI = " << phi[mcp] << std::endl;
490  if (p4v.perp() < 1E-3) phi[mcp] = 0;
491 
492 
493  // just debugging!
494  Hep3Vector gvvec (gvout[0], gvout[1], gvout[2]);
495  float gphi = gvvec.phi(); // gives [-PI, PI]
496  streamlog_out(DEBUG) << "gphi on BeamCal surface in global coordinates = " << gphi << std::endl;
497  if (gphi < 0) gphi += 2*PI; // have [0, 2 PI] now
498  streamlog_out(DEBUG) << "gphi after +2PI = " << gphi << std::endl;
499  float gtheta = gvvec.theta();
500  streamlog_out(DEBUG) << "gtheta on BeamCal surface in global coordinates = " << gtheta << std::endl;
501  streamlog_out(DEBUG) << "radius on BeamCal surface in global coordinates = " << gvvec.perp()*lvvec.mag()/gvvec.mag() << std::endl;
502 
503  // Loop to sum the energy density over all layers.
504  // ebkg = Energy of background
505  // If res == 1 GetEnergyDensity succeeded
506 
507  ebkg[mcp] = 0;
508  ebkg_err[mcp] = 0;
509  bool res = 0;
510  double en_dens = 0;
511  double en_dens_err = 0;
512 
513  for (int layer=1; layer<31; layer++) {
514  if (!newMap && layer ==1) continue; // skip layer if using old map!
515 
516  int sign = 1;
517  if(lvout[2] < 0) sign = -1;
518  res = bc_en->GetEnergyDensity( layer * sign, radius[mcp], phi[mcp], &en_dens, &en_dens_err);
519  streamlog_out(DEBUG) << "energy density in layer " << layer * sign << " = " << en_dens << std::endl;
520 
521  if (res) {
522  ebkg[mcp] += en_dens;
523  ebkg_err[mcp] += en_dens_err*en_dens_err;
524  }
525 
526  }
527  if (ebkg_err[mcp]>=0) {
528  ebkg_err[mcp] = sqrt(ebkg_err[mcp]);
529  }
530  else {
531  ebkg_err[mcp] = 0;
532  }
533  streamlog_out(DEBUG) << "background energy density sum = " << ebkg[mcp] << std::endl;
534 
535 
536  //Scaling of the Energy density with scaling factor
537  ebkg[mcp] *= densityScaling;
538 
539  if (densityScaling != 1.) {
540  streamlog_out(DEBUG) << "scaling background energy density sum with factor " << densityScaling
541  << ", scaled density = " << ebkg[mcp] <<std::endl;
542  }
543 
544 
545  // Here starts the Eff_routine that parametrizes the efficiency
546 
547  //assuming that efficiency is parametrised in GeV/cm^2 instead of GeV/mm^2 ...
548  //alternative: efficiency is in GeV / cell -> *= 55
549  ebkg[mcp] *= 100.;
550 
551  efficiency[mcp] = 0;
552 
553  if ( ebkg[mcp] <= 0 && useInputClusters) { // if == 0, then particle not in BCAL acceptance!
554 
555  streamlog_out(DEBUG) << "Particle made SGV cluster, but ebkg[mcp] = " << ebkg[mcp]
556  << ", so not in map: radius[mcp] = " << radius[mcp]
557  << ", phi[mcp] = " << phi[mcp]*180/3.1415 << std::endl;
558  }
559 
560  else if ( ebkg[mcp] > 0 ) { // if == 0, then particle not in BCAL acceptance!
561 
562  float p0, p1, p2;
563 
564  //energy[mcp] = sqrt(pow(pout[0],2)+pow(pout[1],2)+pow(pout[2],2));
565 
566  if (ebkg[mcp]>=0 && ebkg[mcp]<35 && energy[mcp]>45 && energy[mcp]<345) {
567  if (energy[mcp]>255) {
568  streamlog_out(WARNING) << "not in original range of parametrisation , energy = " << energy[mcp] << std::endl;
569  }
570 
571  p0 = -3.2754e-2 + 8.62e-4*energy[mcp] -2.4424e-6*energy[mcp]*energy[mcp];
572  p1 = -8.936e-3 + 6.05e-4*energy[mcp] -1.719e-6*energy[mcp]*energy[mcp];
573  p2 = 6.112e-2 + 2.84e+3/(energy[mcp]*energy[mcp]);
574 
575  efficiency[mcp] = p0 + 1/(1+p1*exp(p2*ebkg[mcp]));
576 
577  if (efficiency[mcp] > 1.) efficiency[mcp] = 1.;
578  if (efficiency[mcp] < 0.) efficiency[mcp] = 0.;
579 
580  }
581  else if (ebkg[mcp]>=0 && ebkg[mcp]<35 && energy[mcp]>45 && energy[mcp]<550) {
582  efficiency[mcp] = 1.;
583  }
584  else {
585 
586  streamlog_out(ERROR) << "not in valid range of parametrisation , ebkg = "
587  << ebkg[mcp] << ", energy = " << energy[mcp] << std::endl;
588 
589  }
590 
591  // End of the Eff_routine
592 
593  // Selection algorithm
594 
595  char energystring [50];
596  char charseed[5];
597  int seed;
598 
599  // seed setting
600  sprintf(energystring,"%f",energy[mcp]);
601 
602  charseed[0]=energystring[4];
603  charseed[1]=energystring[5];
604  charseed[2]=energystring[6];
605  charseed[3]=energystring[7];
606  charseed[4]=energystring[8];
607 
608  seed = atoi(charseed);
609 
610  TRandom3 *randGen = new TRandom3(seed);
611 
612  rand[mcp] = randGen->Rndm();
613 
614  tag[mcp] = 0;
615  if (rand[mcp] < efficiency[mcp] || (detectAll == 1 && efficiency[mcp]>0.0) ) {
616 
617  ReconstructedParticleImpl* particle = new ReconstructedParticleImpl;
618  ClusterImpl* cluster = new ClusterImpl;
619  streamlog_out(DEBUG) << "adding relation between particle " << particle << "and MCParticle " << p << endl;
620  bcalTruthRelNav.addRelation( particle , p , 1.0) ;
621  //LCRelationImpl* MCrel = new LCRelationImpl(particle,p,0.5);
622  //MCrel->setFrom (particle);
623  //MCrel->setTo (p);
624 
625  //FG: take generator mass ( not necessarily equal to PDG mass )
626  const double m = p->getMass() ;
627  particle->setMass(m);
628  particle->setCharge(q);
629 
630  if (DBDsample) {
631  pxPrime[mcp] = pxIP[mcp]; // no boost needed
632  }
633  else {
634  // DO BOOST!
635  // before the transformation: p4v
636  // (will store 4 vector at IP, not BCAL surface!)
637  // the following is taken from Mokka:
638  // PrimaryGeneratorAction::ApplyLorentzTransformation
639  pxPrime[mcp] = betagamma * energy[mcp] + gamma * pxIP[mcp];
640  }
641  // py and pz remain the same, E changes implicitly with px
642 
643  scaleP[mcp] = 1.;
644 
645  // energy resolution:
646  // dE = 25% * sqrt(E) (intrinsic resolution) + 15% * E_bgk_dens * 55mm^2 (background fluctuations)
647  // alternative: dE = 130% sqrt(E) in total
648  //double sigmaE = 1.3*sqrt(energy[mcp]);
649  //double sigmaE = 0.25*sqrt(energy[mcp]) + 0.15 * ebkg[mcp] * 55; // ????
650  double sigmaE = 0.25*sqrt(energy[mcp]) + ebkg_err[mcp] * 55;
651  if (smearEnergy == true) {
652  double deltaE = randGen->Gaus() * sigmaE;
653  scaleP[mcp] = (1+deltaE/energy[mcp]);
654  if (scaleP[mcp] < 0) scaleP[mcp] = 0;
655  }
656 
657  double pnew[3] = {scaleP[mcp]*pxPrime[mcp], scaleP[mcp]*pyIP[mcp], scaleP[mcp]*pzIP[mcp]};
658  particle->setMomentum(pnew);
659  particle->setEnergy(scaleP[mcp]*sqrt(m*m + pxPrime[mcp]*pxPrime[mcp] +pyIP[mcp]*pyIP[mcp] +pzIP[mcp]*pzIP[mcp]));
660 
661  // store efficiency as "GoodnessOfPID"
662  particle->setGoodnessOfPID(efficiency[mcp]);
663 
664  cluster->setEnergy(particle->getEnergy());
665  cluster->setEnergyError(sigmaE);
666  cluster->setPosition(gvout); // global ILD coord
667  cluster->setITheta(gtheta);
668  cluster->setIPhi(gphi);
669  cluster->subdetectorEnergies().resize(6); cluster->subdetectorEnergies()[5] = particle->getEnergy();
670 
671  particle->addCluster(cluster);
672 
673  BCALCol->addElement(particle);
674  //BCALEffiMCTruthLink->addElement(MCrel);
675  BCALClusters->addElement(cluster);
676 
677  // store positions as single values in tree
678  lposx[mcp]=lvout[0]; // local BeamCal coord
679  lposy[mcp]=lvout[1];
680  lposz[mcp]=lvout[2];
681  gposx[mcp]=gvout[0]; // global ILD coord
682  gposy[mcp]=gvout[1];
683  gposz[mcp]=gvout[2];
684  tag[mcp] = 1;
685 
686 
687  ePrime[mcp] = particle->getEnergy();
688 
689 
690  streamlog_out(DEBUG) << "Detected particle with ID " << pdg[mcp] << endl;
691  streamlog_out(DEBUG) << "Energy before smearing = " << energy[mcp] << endl;
692  streamlog_out(DEBUG) << "Energy after smearing = " << ePrime[mcp] << endl;
693  streamlog_out(DEBUG) << "Efficiency= " << efficiency[mcp] << endl;
694  streamlog_out(DEBUG) << "E background= " << ebkg[mcp] << endl;
695  streamlog_out(DEBUG) << "+++++++++++++++++++" << endl;
696 
697  }
698 
699  mcp++;
700 
701  if (writeTree) {
702  streamlog_out(DEBUG) << "filling tree, mcp = " << mcp << endl;
703  tree->Fill();
704  }
705 
706  delete randGen;
707 
708  } // beamcal acceptance
709 
710  } // energy thresholds
711 
712  } // end of MCparticle loop
713 
714  }
715 
716  BCALEffiMCTruthLink = bcalTruthRelNav.createLCCollection() ;
717  streamlog_out(DEBUG) << "name of BCALEffiMCTruthLink: " << _BCALEffiMCTruthLinkName << ", number of elements = " << BCALEffiMCTruthLink->getNumberOfElements() << endl;
718  streamlog_out(DEBUG) << "name of BCALCol: " << _BCALParticleName << ", number of elements = " << BCALCol->getNumberOfElements() << endl;
719 
720  //if ( BCALCol->getNumberOfElements() <= 0 ) {
721  //delete BCALCol;
722  //delete BCALEffiMCTruthLink;
723  //delete BCALClusters;
724  //} else {
725  evt->addCollection(BCALCol ,_BCALParticleName) ;
726  evt->addCollection(BCALEffiMCTruthLink ,_BCALEffiMCTruthLinkName) ;
727  evt->addCollection(BCALClusters ,_BCALClusterName) ;
728  //}
729 
730 
731  streamlog_out(DEBUG) << "+++++++++++++++++ finished processing event: " << evt->getEventNumber()
732  << " in run: " << evt->getRunNumber()
733  << std::endl ;
734 
735 
736  _nEvt ++ ;
737 }
738 
739 
740 
741 void BCalTagEfficiency::check( LCEvent * evt ) {
742  // nothing to check here - could be used to fill checkplots in reconstruction processor
743 }
744 
745 
746 void BCalTagEfficiency::end(){
747 
748  delete bc_en;
749 
750  if (writeTree) {
751  tree->Write();
752  rootfile->Write();
753  }
754 
755 }
756 
BCalTagEfficiency::_BCALInputTruthLinkName
std::string _BCALInputTruthLinkName
Definition: BCalTagEfficiency.h:81
BcEnergyDensity::GetEnergyDensity
Bool_t GetEnergyDensity(const Int_t &rLayer, const Double_t &rRadius, const Double_t &rPhi, Double_t *pEnDens, Double_t *pEnDensError) const
Definition: BcEnergyDensity.cc:99
BCalTagEfficiency::end
virtual void end()
Called after data processing for clean up.
Definition: BCalTagEfficiency.cc:746
BCalTagEfficiency::processEvent
virtual void processEvent(LCEvent *evt)
Called for every event - the working horse.
Definition: BCalTagEfficiency.cc:327
aBCalTagEfficiency
BCalTagEfficiency aBCalTagEfficiency
Definition: BCalTagEfficiency.cc:58
BCalTagEfficiency::backgroundfilename
std::string backgroundfilename
Definition: BCalTagEfficiency.h:111
sistrip::m
static const float m
Definition: PhysicalConstants.h:36
Vector3D
CLHEP::Hep3Vector Vector3D
Definition: DistilledPFOCreator.cc:11
BCalTagEfficiency::_BCALEffiMCTruthLinkName
std::string _BCALEffiMCTruthLinkName
Definition: BCalTagEfficiency.h:82
BCalTagEfficiency::_BCALClusterName
std::string _BCALClusterName
Definition: BCalTagEfficiency.h:83
BCalTagEfficiency::check
virtual void check(LCEvent *evt)
Definition: BCalTagEfficiency.cc:741
BCalTagEfficiency
Marlin processor to calculate BCAL tagging efficiency.
Definition: BCalTagEfficiency.h:45
BCalTagEfficiency::_MCParticleName
std::string _MCParticleName
Input collection name.
Definition: BCalTagEfficiency.h:80
BCalTagEfficiency::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Called for every run.
Definition: BCalTagEfficiency.cc:322
sistrip::LCCollectionVec
std::vector< LCCollection * > LCCollectionVec
Definition: SiStripClus.h:55
BCalTagEfficiency::_BCALParticleName
std::string _BCALParticleName
Definition: BCalTagEfficiency.h:84
BCalTagEfficiency::init
virtual void init()
Called at the begin of the job before anything is read.
Definition: BCalTagEfficiency.cc:179
BCalTagEfficiency::bc_en
BcEnergyDensity * bc_en
Definition: BCalTagEfficiency.h:130
BCalTagEfficiency::SGVmapfilename
std::string SGVmapfilename
Definition: BCalTagEfficiency.h:112
bcalhit_
void bcalhit_(float pin_[3], float &q_, float vin_[3], float &B_, float &eBeam_, float &zbcal_, float pout_[3], float vout_[3])