All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
AbsCalibr.cc
Go to the documentation of this file.
1 #include "AbsCalibr.h"
2 #include "Calibration.h"
3 
4 #include <cstdio>
5 #include <cmath>
6 #include <cstring>
7 #include "lcio.h"
8 #include <marlin/Global.h>
9 #include <gear/GEAR.h>
10 #include <gear/CalorimeterParameters.h>
11 #include "EVENT/LCCollection.h"
12 #include "EVENT/SimCalorimeterHit.h"
13 #include "EVENT/CalorimeterHit.h"
14 #include "EVENT/RawCalorimeterHit.h"
15 #include "EVENT/MCParticle.h"
16 #include "IMPL/LCCollectionVec.h"
17 
18 
19 #include <cstdlib>
20 
21 // To create collection can be output by LCIO
22 #include "UTIL/LCFixedObject.h"
23 
24 // #define ASCII_OUTPUT
25 // #define USE_HBOOK
26 
27 #ifdef USE_HBOOK
28 //+++++++++++++ HBOOK +++ HPLOT ++++++++++++++++
29 #ifndef f2cFortran
30 #define f2cFortran 1
31 #endif
32 #include <cfortran.h>
33 #include <graflib.h>
34 //#include <packlib.h>
35 //#include <kernlib.h>
36 #include <hbook.h>
37 #include <higz.h>
38 #define PAWC_SIZE 5000000
39 typedef struct { float PAW[PAWC_SIZE]; } PAWC_DEF;
40 #define PAWC COMMON_BLOCK(PAWC,pawc)
41 COMMON_BLOCK_DEF(PAWC_DEF,PAWC);
42 //+++++++++++++ HBOOK +++ HPLOT ++++++++++++++++
43 #endif
44 
45 
46 // using namespace std ;
47 using namespace EVENT ;
48 using namespace IMPL ;
49 using namespace UTIL;
50 
51 #ifdef USE_HBOOK
52 // +++++++++++++ HBOOK +++ HPLOT ++++++++++++++++
53 //============================================================================
54 static void plot_init(){ // Open HIGZ windows and Initialize HBOOK
55 //============================================================================
56  static bool is_init=0;
57  if(is_init) return;
58  HLIMIT(PAWC_SIZE);
59  HPLINT(4);
60  is_init=1;
61 } // End of plot_init()
62 
63 //============================================================================
64 static void hist_init(){ // Book Histograms
65 //============================================================================
66  static bool is_init=0;
67  if(is_init) return;
68 
69 
70  // 1000 GeV ECM
71 // HBOOK1(1,"ECAL hit energy sum ",700,0.0,1200.,0.);
72 // HBOOK1(2,"HCAL hit energy sum ",700,0.0,1200.,0.);
73 // HBOOK1(3,"Calorimeter hit energy sum ",700,0.0,1200.,0.);
74 // HBOOK2(4,"E-ECAL vs E-HCAL",125,0.0,1000.0, 125,0.0,1000.0, 0.0);
75 // HBOOK1(5,"Calorimeter hit energy sum ",350,0.0,1200.,0.);
76 // HBOOK1(6,"Calorimeter energy - Available",120,-120.0,120.,0.);
77 //!!!
78  // 500 GeV ECM
79  HBOOK1(1,"ECAL hit energy sum ",700,0.0,700.,0.);
80  HBOOK1(2,"HCAL hit energy sum ",700,0.0,700.,0.);
81  HBOOK1(3,"Calorimeter hit energy sum ",700,0.0,700.,0.);
82  HBOOK2(4,"E-ECAL vs E-HCAL",125,0.0,500.0, 125,0.0,500.0, 0.0);
83  HBOOK1(5,"Calorimeter hit energy sum ",350,0.0,700.,0.);
84  HBOOK1(6,"Calorimeter energy - Available",120,-60.0,60.,0.);
85 
86  // HPLOPT("GRID",1);
87  // HPLOPT("STAT",1);
88  is_init=1;
89 } //End of hist_init()
90 // +++++++++++++ HBOOK +++ HPLOT ++++++++++++++++
91 #endif
92 
93 namespace marlin{
94 
95 double Balance(LCEvent * evt);
96 
97  AbsCalibr aAbsCalibr ;
98  AbsCalibr::AbsCalibr() : Processor("AbsCalibr") {
99  _description = " Calorimeter Abslute Energy Calibration" ;
100  vector<int> nlayer;
101  nlayer.push_back(30);
102  nlayer.push_back(10);
103  nlayer.push_back(40);
104  registerProcessorParameter("NLayer","Number of layers in zone",_nlayer,nlayer);
105  vector<float> coeff;
106  coeff.push_back(33.02346);
107  coeff.push_back(93.56822);
108  coeff.push_back(21.196262);
109  registerProcessorParameter("Coeff","Calorimeter coeffs",_coeff,coeff);
110  vector<float> cuts;
111  cuts.push_back(170.0e-6 * 0.5);
112  cuts.push_back(170.0e-6 * 0.5);
113  cuts.push_back(875.0e-6 * 0.5);
114  registerProcessorParameter("Cuts","Calorimeter cuts",_cuts,cuts);
115  }
116 
117 //============================================================================
118  void AbsCalibr::init() {
119 //============================================================================
120  std::cout << "AbsCalibr::init() " << name() << std::endl
121  << " parameters: " << std::endl << *parameters() ;
122 
123 
124 
125 #ifdef USE_HBOOK
126 // +++++++++++++ HBOOK +++ HPLOT ++++++++++++++++
127  plot_init();
128  hist_init();
129 
130  std::cout <<
131  "++++++++++++++ HBOOK and HPLOT were initiaslized +++++++++++++++++++ "
132  << std::endl;
133 // +++++++++++++ HBOOK +++ HPLOT ++++++++++++++++
134 #endif
135 
136  _nRun = 0 ;
137  _nEvt = 0 ;
138 }
139 
140 //============================================================================
141  void AbsCalibr::processRunHeader( LCRunHeader* run) {
142 //============================================================================
143  std::cout << "AbsCalibr::processRun() " << name()
144  << " in run " << run->getRunNumber() << "; description is "
145  <<run->getDescription()<<std::endl ;
146  _nRun++ ;
147  }
148 //============================================================================
149  void AbsCalibr::processEvent( LCEvent * evt ) {
150 //============================================================================
151 
152  LCCollectionVec* abc = new LCCollectionVec(LCIO::LCGENERICOBJECT);
153 
154  CalorimeterHit* hit; // pointer of class CalorimeterHit
155 
156  int nHits[3];
157  double en[3];
158 
159  nHits[ECAL1]=nHits[ECAL2]=nHits[HCAL]=0;
160  en[ECAL1]=en[ECAL2]=en[HCAL]=0.;
161 
162 // Find collection names to fill hits arrayes
163  const std::vector<std::string>* strVec = evt->getCollectionNames() ;
164  for( std::vector<std::string>::const_iterator name = strVec->begin() ;
165  name != strVec->end() ; name++) {
166  const std::string & tname = evt->getCollection( *name )->getTypeName();
167  // std::cout<<"Collection type: "<<tname<<std::endl;
168 
169 //++++++++++++++++++++++ ECAL hit filling +++++++++++++++++++++++++++++++++++++++
170  if(tname == LCIO::CALORIMETERHIT ){
171 //-----------------------------------------------------------------------
172  if(!std::strncmp((*name).data(),"ECAL",4)){ // Resolve collection name
173 //-----------------------------------------------------------------------
174 // std::cout << " ECAL: " << *name << std::endl;
175  std::basic_string<char> coll = *name;
176  EVENT::LCCollection* col_ECAL = evt->getCollection(coll) ;
177 
178  unsigned n = col_ECAL->getNumberOfElements(); // Get number of hits in ECAL
179 
180  for( unsigned i = 0 ; i < n ; i++ ){ // then fill it
181  hit = dynamic_cast<CalorimeterHit*>( col_ECAL->getElementAt(i));
182  float e = hit->getEnergy();
183  unsigned keyGE = hit->getCellID0() ;
184  int l = (keyGE & 0x3F000000)>>24;
185  if(l<_nlayer[ECAL1]){
186  if(e>_cuts[ECAL1]){
187  en[ECAL1]+=e*_coeff[ECAL1];
188  nHits[ECAL1]++;
189  }
190  } else {
191  if(e>_cuts[ECAL2]){
192  en[ECAL2]+=e*_coeff[ECAL2];
193  nHits[ECAL2]++;
194  }
195  } // end ECAL hit loop
196  } // Resolve collection name
197 
198 //+++++++++++++++++++++ HCAL hit filling +++++++++++++++++++++++++++++++++++++++++++++
199  } else if(!std::strncmp((*name).data(),"HCAL",4)){// Resolve collection name
200 //-----------------------------------------------------------------------
201 // std::cout << " HCAL: " << *name << std::endl;
202  std::basic_string<char> coll = *name;
203  EVENT::LCCollection* col_HCAL = evt->getCollection(coll) ;
204 
205  nHits[HCAL] = col_HCAL->getNumberOfElements() ; //Get number of hits in HCAL
206 
207  for( int i=0 ; i< nHits[HCAL] ; i++ ){ // then fill it
208  hit = dynamic_cast<CalorimeterHit*>( col_HCAL->getElementAt(i));
209  float e = hit->getEnergy();
210  if(e>_cuts[HCAL])
211  en[HCAL]+=e*_coeff[HCAL];
212  } // end HCAL hit loop
213  } // Resolve collection name
214  else
215  std::cout << " !!! UNKNOWN !!!: " << *name << std::endl;
216  } // along calorimeter collections in LCIO
217  } // along all collections in LCIO
218 
219  double E_Real = Balance(evt);
220 
221 //-----------------------------------------------------------------------
222 // Create collection to store it in *.slcio
223 //-----------------------------------------------------------------------
224  Calibration* b = new Calibration(_nEvt,nHits[0],nHits[1],nHits[2],en[0],en[1],en[2],E_Real);
225  abc->addElement(b);
226  evt->addCollection(abc, "AbsCalibr");
227 //-----------------------------------------------------------------------
228 
229 #ifdef ASCII_OUTPUT
230  static FILE *f=NULL;
231  if(!f) f=fopen("abs_calibr.root","w");
232  fprintf(f,"%10i %10u %10u %10u %15.3f %15.3f %15.3f %15.3f\n",
233  _nEvt,nHits[ECAL1],nHits[ECAL2],nHits[HCAL],
234  en[ECAL1],en[ECAL2],en[HCAL],E_Real);
235 #endif
236 
237 
238 
239 #ifdef USE_HBOOK
240 /*******
241  if (E_Real> 0.5) {
242  HF1(1,E_Ecal,1.0);
243  HF1(2,E_Hcal,1.0);
244  HF1(3,E_Ecal+E_Hcal,1.0);
245  HF1(5,E_Ecal+E_Hcal,1.0);
246  HF2(4,E_Ecal,E_Hcal,1.0);
247  HF1(6,E_Ecal+E_Hcal-E_Real,1.0);
248  }
249 
250  if((_nEvt%20) == 0){
251  // Draw current Histograms
252  HPLOPT("STAT",1);
253  // HPLOPT("LOGY",1);
254  HPLZON(2,2,1," ");
255  HPLOT(5,"BOX ","HIST",0);
256  HPLOT(6,"BOX ","HIST",0);
257  HPLOT(3,"BOX ","HIST",0);
258  HPLOT(4,"BOX ","HIST",0);
259  IUWK(0,1);
260  HPLOPT("NSTA",1);
261  }
262 ******/
263 #endif
264  // std::cout<< " ---> Calo energy = " << E_Ecal <<" + "<< E_Hcal
265  // <<" = "<< E_Ecal+E_Hcal<<"; Available E = "<<E_Real<<std::endl;
266 
267 #ifdef USE_HBOOK
268 // Wait for next event
269 //std::cout << " ++++++++++ TYPE o TO SAVE HISTOGRAMS, IF YOU WISH, +++++++" << std::endl ;
270  // int c = getchar();
271  // if(c=='o') HRPUT(0,"abs_calibr.his"," "); // Save all existing histograms on disk
272 #endif
273 
274  _nEvt ++ ;
275  return ; // from void AbsCalibr::processEvent( LCEvent * evt )
276 }
277 
278 //============================================================================
279 void AbsCalibr::check( LCEvent * /*evt*/ ) { ;}
280 
281 //============================================================================
282 void AbsCalibr::end(){
283 //============================================================================
284  std::cout << "AbsCalibr::end() " << name()
285  << " processed " << _nEvt << " events in " << _nRun << " runs "
286  << std::endl ;
287 #ifdef USE_HBOOK
288  HRPUT(0,"abs_calibr.his"," "); // Save all existing histograms on disk
289 #endif
290 }
291 
292 //============================================================================
293 double Balance(LCEvent * evt){
294 //============================================================================
295  int idpdg;
296  const double* mom;
297  float enr;
298  double mass;
299  LCCollection* mcpCol = evt->getCollection("MCParticle" ) ;
300 //-----------------------------------------------------------------------
301 // Calculate balance at IP taking into account everything
302 //-----------------------------------------------------------------------
303  double px,py,pz,pt,ttet;
304 
305  double e_to_tube = 0.;
306  double e_to_tubex = 0.;
307  double e_to_tubey = 0.;
308  double e_to_tubez = 0.;
309  int n_to_tube = 0;
310 
311  double e_neutr = 0.;
312  double e_neutrx= 0.;
313  double e_neutry= 0.;
314  double e_neutrz= 0.;
315  int n_neutr= 0;
316 
317  double e_muon = 0.;
318  double e_muonx= 0.;
319  double e_muony= 0.;
320  double e_muonz= 0.;
321  int n_muon= 0;
322 
323  double e_elect = 0.;
324  double e_electx= 0.;
325  double e_electy= 0.;
326  double e_electz= 0.;
327  int n_elect= 0;
328 
329  double e_photon = 0.;
330  double e_photonx= 0.;
331  double e_photony= 0.;
332  double e_photonz= 0.;
333  int n_photon= 0;
334 
335  double e_pi0 = 0.;
336  double e_pi0x= 0.;
337  double e_pi0y= 0.;
338  double e_pi0z= 0.;
339  int n_pi0= 0;
340 
341  double e_llhadr = 0.;
342  double e_llhadrx= 0.;
343  double e_llhadry= 0.;
344  double e_llhadrz= 0.;
345  int n_llhadr= 0;
346 
347  double e_slhadr = 0.;
348  double e_slhadrx= 0.;
349  double e_slhadry= 0.;
350  double e_slhadrz= 0.;
351  int n_slhadr= 0;
352 
353  double e_chadr = 0.;
354  double e_chadrx= 0.;
355  double e_chadry= 0.;
356  double e_chadrz= 0.;
357  int n_chadr= 0;
358 
359  for(int i=0; i<mcpCol->getNumberOfElements() ; i++){
360  MCParticle* imc = dynamic_cast<MCParticle*> ( mcpCol->getElementAt( i ) ) ;
361  idpdg = imc-> getPDG ();
362  mom = imc-> getMomentum ();
363  enr = imc-> getEnergy ();
364  mass = imc-> getMass ();
365  if( imc-> getGeneratorStatus() == 1) { // stable particles only
366  px = mom[0];
367  py = mom[1];
368  pz = mom[2];
369  pt = hypot(px,py);
370  ttet = atan2(pt,pz);
371  if ((fabs(ttet) < 0.1) || (fabs(M_PI-ttet) < 0.1)) {
372  e_to_tube += enr;
373  e_to_tubex += px;
374  e_to_tubey += py;
375  e_to_tubez += pz;
376  n_to_tube ++;
377  continue;
378  }
379  if((abs(idpdg)==12)||(abs(idpdg)==14)||(abs(idpdg)==16)) {
380  e_neutr += enr;
381  e_neutrx += px;
382  e_neutry += py;
383  e_neutrz += pz;
384  n_neutr ++;
385  continue;
386  }
387  if(abs(idpdg)==13) { // mu+ mu-
388  e_muon += enr;
389  e_muonx += px;
390  e_muony += py;
391  e_muonz += pz;
392  n_muon ++;
393  continue;
394  }
395  if(abs(idpdg)==11) { // e+ e-
396  e_elect += enr;
397  e_electx += px;
398  e_electy += py;
399  e_electz += pz;
400  n_elect ++;
401  continue;
402  }
403  if(idpdg == 111) { // Pi0 as stable
404  e_pi0 += enr;
405  e_pi0x += px;
406  e_pi0y += py;
407  e_pi0z += pz;
408  n_pi0 ++;
409  continue;
410  }
411  if(idpdg == 22) { // photon
412  e_photon += enr;
413  e_photonx += px;
414  e_photony += py;
415  e_photonz += pz;
416  n_photon ++;
417  continue;
418  }
419  if( // long lived neutral hadrons
420  (abs(idpdg)==2112)|| // neutron
421  (abs(idpdg)== 130) // KoL
422  ) {
423  e_llhadr += enr;
424  e_llhadrx += px;
425  e_llhadry += py;
426  e_llhadrz += pz;
427  n_llhadr ++;
428  continue;
429  }
430  if( // short lived neutral hadrons
431  (abs(idpdg)== 310)|| // KoS
432  (abs(idpdg)==3122)|| // Lambda0
433  (abs(idpdg)==3212)|| // Sigma0
434  (abs(idpdg)==3322) // Xi0
435  ) {
436  e_slhadr += enr;
437  e_slhadrx += px;
438  e_slhadry += py;
439  e_slhadrz += pz;
440  n_slhadr ++;
441  continue;
442  }
443  if(!(abs(idpdg)==12) && !(abs(idpdg)==14) && !(abs(idpdg)==16) && // neutrinos
444  !(abs(idpdg)==13) && // mu+ mu-
445  !(abs(idpdg)==11) && // e+ e-
446  !(idpdg == 111) && // Pi0
447  !(idpdg == 22) && // photon
448  !(abs(idpdg)==2112) && !(abs(idpdg)== 311) && // neutral hadrons
449  !(abs(idpdg)== 130) && !(abs(idpdg)== 310) && // neutral hadrons
450  !(abs(idpdg)==3122) && !(abs(idpdg)==3212)) { // neutral hadrons
451  e_chadr += enr;
452  e_chadrx += px;
453  e_chadry += py;
454  e_chadrz += pz;
455  n_chadr ++;
456  continue;
457  }
458  std::cout <<" Unknow for this program ID is " <<idpdg<< std::endl;
459  } // Stable particles only
460  } // End for for MCParticles
461  double e_sum = e_elect+e_muon+e_chadr+e_pi0+e_photon+e_llhadr+e_slhadr+e_neutr;
462  double evt_energy = e_sum;
463  // Minus muon loosed energy = 1.3 GeV in average
464  double e_mu_lost = e_muon - n_muon*1.3;
465  double e_lost = e_neutr + e_mu_lost + e_to_tube;
466  double e_real = evt_energy - e_lost;
467  if(e_real< 0.0) e_real = 0.000001;
468 
469  /**
470  std::cout <<" =============================================================="<< std::endl;
471  std::cout << " ======== Record Balance ======="<< std::endl ;
472  std::cout <<" =============================================================="<< std::endl;
473  std::cout <<" ============== Possible lost ==================="<< std::endl;
474  std::cout <<" Neutrino energy = "<<e_neutr<<", in "<< n_neutr<<" neutrinos"<< std::endl;
475  std::cout <<" Energy to beam tube = "<<e_to_tube<<", in "<<n_to_tube<<" particles"<< std::endl;
476  std::cout <<" Muons energy lost = "<<e_mu_lost<<" in "<<n_muon<<" muons"<< std::endl;
477  std::cout <<" --------------------------------------------------"<< std::endl;
478  std::cout <<" Total Event energy at IP = "<<evt_energy<<" [GeV]"<< std::endl;
479  std::cout <<" --------------------------------------------------"<< std::endl;
480  std::cout <<" Whole lost Energy = "<<e_lost<< std::endl;
481  std::cout <<" Available Energy in calo.= "<<e_real<< std::endl;
482  std::cout <<" =============================================================="<< std::endl;
483  std::cout <<" Muon energy = "<<e_muon <<'\t'<<" in "<< n_muon <<" muons"<< std::endl;
484  std::cout <<" Electron energy = "<<e_elect <<'\t'<<" in "<< n_elect <<" electrons"<< std::endl;
485  std::cout <<" Charged hadron energy = "<<e_chadr <<'\t'<<" in "<< n_chadr <<" hadrons"<< std::endl;
486  std::cout <<" -------------------------------------------------------------"<< std::endl;
487  std::cout <<" Pi0 energy (if stable) = "<<e_pi0 <<'\t'<<" in "<< n_pi0 <<" Pi zeros"<< std::endl;
488  std::cout <<" Photon energy = "<<e_photon<<'\t'<<" in "<< n_photon<<" photons"<< std::endl;
489  std::cout <<" -------------------------------------------------------------"<< std::endl;
490  std::cout <<" Long lived hadron energy = "<<e_llhadr<<'\t'<<" in "<< n_llhadr<<" hadrons"<< std::endl;
491  std::cout <<" Short lived hadron energy = "<<e_slhadr<<'\t'<<" in "<< n_slhadr<<" hadrons"<< std::endl;
492  std::cout <<" =============================================================="<< std::endl;
493  */
494 
495  return e_real;
496 
497  } // End MC_Balance
498 
499 }// namespace marlin
marlin::Balance
double Balance(LCEvent *evt)
Definition: AbsCalibr.cc:293
marlin::AbsCalibr::init
virtual void init()
Definition: AbsCalibr.cc:118
marlin::AbsCalibr::ECAL1
Definition: AbsCalibr.h:52
marlin::AbsCalibr::_coeff
vector< float > _coeff
Definition: AbsCalibr.h:58
marlin::AbsCalibr::processEvent
virtual void processEvent(LCEvent *evt)
Definition: AbsCalibr.cc:149
TNT::hypot
Real hypot(const Real &a, const Real &b)
Definition: tnt_math_utils.h:18
marlin::AbsCalibr::_nRun
int _nRun
Definition: AbsCalibr.h:48
marlin::aAbsCalibr
AbsCalibr aAbsCalibr
Definition: AbsCalibr.cc:97
marlin::AbsCalibr::_nEvt
int _nEvt
Definition: AbsCalibr.h:49
sistrip::e
static const float e
Definition: PhysicalConstants.h:26
sistrip::LCCollectionVec
std::vector< LCCollection * > LCCollectionVec
Definition: SiStripClus.h:55
marlin::AbsCalibr::processRunHeader
virtual void processRunHeader(LCRunHeader *run)
Definition: AbsCalibr.cc:141
marlin::AbsCalibr::HCAL
Definition: AbsCalibr.h:54
marlin::AbsCalibr::ECAL2
Definition: AbsCalibr.h:53
marlin::AbsCalibr::_cuts
vector< float > _cuts
Definition: AbsCalibr.h:59
marlin::AbsCalibr::end
virtual void end()
Definition: AbsCalibr.cc:282
marlin::AbsCalibr::_nlayer
vector< int > _nlayer
Definition: AbsCalibr.h:57
marlin::AbsCalibr::check
virtual void check(LCEvent *evt)
Definition: AbsCalibr.cc:279