HCal-analysis/HCALAnalysis/HCALAnalysis.C

0001 //  General sPHENIX tools
0002
0003 #include "HCALAnalysis.h"
0004
0005 #include <calobase/RawTowerContainer.h>
0006 #include <pdbcalbase/PdbParameterMap.h>
0007 #include <phparameter/PHParameters.h>
0008 #include <ffaobjects/EventHeader.h>
0009 #include <g4detectors/PHG4ScintillatorSlat.h>
0010 #include <g4detectors/PHG4ScintillatorSlatContainer.h>
0011 #include <g4detectors/PHG4ScintillatorSlatDefs.h>
0012
0013 #include <fun4all/SubsysReco.h>
0014 #include <fun4all/Fun4AllServer.h>
0015 #include <fun4all/PHTFileServer.h>
0016 #include <phool/PHCompositeNode.h>
0017 #include <fun4all/Fun4AllReturnCodes.h>
0018 #include <phool/getClass.h>
0019 #include <fun4all/Fun4AllHistoManager.h>
0020 #include <phool/getClass.h>
0021
0022
0023 #include <g4main/PHG4HitContainer.h>
0024 #include <g4main/PHG4Hit.h>
0025
0026 //  Root classes
0027 #include <TH1.h>
0028 #include <TH2.h>
0029 #include <TPythia6.h>
0030 #include <TMath.h>
0031 #include <TFile.h>
0032 #include <TNtuple.h>
0033 #include <TLorentzVector.h>
0034
0035 #include <assert.h>
0036
0037
0038 HCALAnalysis::HCALAnalysis(const char *outfile) : SubsysReco("HCAL Analysis Development Code")
0039 {
0040   OutFileName = outfile;
0041
0042   return;
0043 }
0044
0045 int HCALAnalysis::Init(PHCompositeNode *topNode)
0046 {
0047
0048   outputfile = new TFile(OutFileName,"RECREATE");
0049
0050   calenergy = new TNtuple("calenergy","calenergy","e_hcin:e_hcout:e_cemc:ea_hcin:ea_hcout:ea_cemc:ev_hcin:ev_hcout:ev_cemc:e_magnet:e_bh:e_emcelect:e_hcalin_spt:sfemc:sfihcal:sfohcal:sfvemc:sfvihcal:sfvohcal:LCG:emc_LCG:hcalin_LCG:hcalout_LCG:e_svtx:e_svtx_support:e_bh_plus:e_bh_minus");
0051
0052   return 0;
0053 }
0054
0055 int HCALAnalysis::process_event(PHCompositeNode *topNode)
0056 {
0057   GetNodes(topNode);
0058
0059   /*
0060   std::cout << " This event has " << _hcalout_hit_container->size() << " HCALOUT G4 hits" << std::endl;
0061   std::cout << " This event has " << _hcalin_hit_container->size() << " HCALIN G4 hits" << std::endl;
0062   std::cout << " This event has " << _cemc_hit_container->size() << " CEMC G4 hits" << std::endl;
0063
0064   std::cout << " This event has " << _hcalout_abs_hit_container->size() << " ABSORBER_HCALOUT G4 hits" << std::endl;
0065   std::cout << " This event has " << _hcalin_abs_hit_container->size() << " ABSORBER_HCALIN G4 hits" << std::endl;
0066   std::cout << " This event has " << _cemc_abs_hit_container->size() << " ABSORBER_CEMC G4 hits" << std::endl;
0067   */
0068
0069   // Sum up energy in scintillator, absorber, and light yield for
0070   // the three calorimeter sections.
0071
0072   float e_hcin = 0.0, e_hcout = 0.0, e_cemc = 0.0;
0073   float ea_hcin = 0.0, ea_hcout = 0.0, ea_cemc = 0.0;
0074   float ev_hcin = 0.0, ev_hcout = 0.0, ev_cemc = 0.0;
0075   float e_magnet = 0.0, e_bh = 0.0;
0076   float e_emcelect = 0.0, e_hcalin_spt = 0.0;
0077   float e_svtx = 0.0, e_svtx_support = 0.0;
0078   float e_bh_plus = 0.0, e_bh_minus = 0.0;
0079
0080   float lcg = 0.0; //longitudinal center of gravity
0081   float emc_lcg = 0.0;
0082   float hcalin_lcg = 0.0;
0083   float hcalout_lcg = 0.0;
0084
0085   PHG4HitContainer::ConstRange hcalout_hit_range = _hcalout_hit_container->getHits();
0086   for (PHG4HitContainer::ConstIterator hit_iter = hcalout_hit_range.first;
0087        hit_iter != hcalout_hit_range.second; hit_iter++){
0088
0089     PHG4Hit *this_hit =  hit_iter->second;
0090     assert(this_hit);
0091
0092     e_hcout += this_hit->get_edep();
0093     ev_hcout += this_hit->get_light_yield();
0094
0095     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0096     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0097
0098     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0099     hcalout_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0100
0101
0102   }
0103
0104   PHG4HitContainer::ConstRange hcalin_hit_range = _hcalin_hit_container->getHits();
0105   for (PHG4HitContainer::ConstIterator hit_iter = hcalin_hit_range.first;
0106        hit_iter != hcalin_hit_range.second; hit_iter++){
0107
0108     PHG4Hit *this_hit =  hit_iter->second;
0109     assert(this_hit);
0110
0111     e_hcin += this_hit->get_edep();
0112     ev_hcin += this_hit->get_light_yield();
0113
0114     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0115     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0116
0117     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0118     hcalin_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0119
0120   }
0121
0122   PHG4HitContainer::ConstRange cemc_hit_range = _cemc_hit_container->getHits();
0123   for (PHG4HitContainer::ConstIterator hit_iter = cemc_hit_range.first;
0124        hit_iter != cemc_hit_range.second; hit_iter++){
0125
0126     PHG4Hit *this_hit =  hit_iter->second;
0127     assert(this_hit);
0128
0129     e_cemc += this_hit->get_edep();
0130     ev_cemc += this_hit->get_light_yield();
0131
0132     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0133     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0134
0135     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0136     emc_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0137
0138   }
0139
0140   PHG4HitContainer::ConstRange hcalout_abs_hit_range = _hcalout_abs_hit_container->getHits();
0141   for (PHG4HitContainer::ConstIterator hit_iter = hcalout_abs_hit_range.first;
0142        hit_iter != hcalout_abs_hit_range.second; hit_iter++){
0143
0144     PHG4Hit *this_hit =  hit_iter->second;
0145     assert(this_hit);
0146
0147     ea_hcout += this_hit->get_edep();
0148
0149     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0150     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0151
0152     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0153     hcalout_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0154
0155   }
0156
0157   PHG4HitContainer::ConstRange hcalin_abs_hit_range = _hcalin_abs_hit_container->getHits();
0158   for (PHG4HitContainer::ConstIterator hit_iter = hcalin_abs_hit_range.first;
0159        hit_iter != hcalin_abs_hit_range.second; hit_iter++){
0160
0161     PHG4Hit *this_hit =  hit_iter->second;
0162     assert(this_hit);
0163
0164     ea_hcin += this_hit->get_edep();
0165
0166     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0167     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0168
0169     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0170     hcalin_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0171
0172   }
0173
0174   PHG4HitContainer::ConstRange cemc_abs_hit_range = _cemc_abs_hit_container->getHits();
0175   for (PHG4HitContainer::ConstIterator hit_iter = cemc_abs_hit_range.first;
0176        hit_iter != cemc_abs_hit_range.second; hit_iter++){
0177
0178     PHG4Hit *this_hit =  hit_iter->second;
0179     assert(this_hit);
0180
0181     ea_cemc += this_hit->get_edep();
0182
0183     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0184     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0185
0186     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0187     emc_lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0188
0189   }
0190
0191   PHG4HitContainer::ConstRange magnet_hit_range = _magnet_hit_container->getHits();
0192   for (PHG4HitContainer::ConstIterator hit_iter = magnet_hit_range.first;
0193        hit_iter != magnet_hit_range.second; hit_iter++){
0194
0195     PHG4Hit *this_hit =  hit_iter->second;
0196     assert(this_hit);
0197
0198     e_magnet += this_hit->get_edep();
0199
0200     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0201     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0202
0203     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0204
0205   }
0206
0207   PHG4HitContainer::ConstRange bh_hit_range = _bh_hit_container->getHits();
0208   for (PHG4HitContainer::ConstIterator hit_iter = bh_hit_range.first;
0209        hit_iter != bh_hit_range.second; hit_iter++){
0210
0211     PHG4Hit *this_hit =  hit_iter->second;
0212     assert(this_hit);
0213
0214     e_bh += this_hit->get_edep();
0215
0216   }
0217
0218   PHG4HitContainer::ConstRange bh_plus_hit_range = _bh_plus_hit_container->getHits();
0219   for (PHG4HitContainer::ConstIterator hit_iter = bh_plus_hit_range.first;
0220        hit_iter != bh_plus_hit_range.second; hit_iter++){
0221
0222     PHG4Hit *this_hit =  hit_iter->second;
0223     assert(this_hit);
0224
0225     e_bh_plus += this_hit->get_edep();
0226
0227   }
0228
0229   PHG4HitContainer::ConstRange bh_minus_hit_range = _bh_minus_hit_container->getHits();
0230   for (PHG4HitContainer::ConstIterator hit_iter = bh_minus_hit_range.first;
0231        hit_iter != bh_minus_hit_range.second; hit_iter++){
0232
0233     PHG4Hit *this_hit =  hit_iter->second;
0234     assert(this_hit);
0235
0236     e_bh_minus += this_hit->get_edep();
0237
0238   }
0239
0240   PHG4HitContainer::ConstRange cemc_electronics_hit_range = _cemc_electronics_hit_container->getHits();
0241   for (PHG4HitContainer::ConstIterator hit_iter = cemc_electronics_hit_range.first;
0242        hit_iter != cemc_electronics_hit_range.second; hit_iter++){
0243
0244     PHG4Hit *this_hit =  hit_iter->second;
0245     assert(this_hit);
0246
0247     e_emcelect += this_hit->get_edep();
0248
0249     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0250     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0251
0252     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0253
0254   }
0255
0256   PHG4HitContainer::ConstRange hcalin_spt_hit_range = _hcalin_spt_hit_container->getHits();
0257   for (PHG4HitContainer::ConstIterator hit_iter = hcalin_spt_hit_range.first;
0258        hit_iter != hcalin_spt_hit_range.second; hit_iter++){
0259
0260     PHG4Hit *this_hit =  hit_iter->second;
0261     assert(this_hit);
0262
0263     e_hcalin_spt += this_hit->get_edep();
0264
0265     float rin =  sqrt(pow(this_hit->get_x(0),2)+pow(this_hit->get_y(0),2));
0266     float rout =  sqrt(pow(this_hit->get_x(1),2)+pow(this_hit->get_y(1),2));
0267
0268     lcg += this_hit->get_edep()*(rin + ((rout-rin)/2.0));
0269
0270   }
0271
0272   // SVTX
0273
0274   PHG4HitContainer::ConstRange svtx_hit_range = _svtx_hit_container->getHits();
0275   for (PHG4HitContainer::ConstIterator hit_iter = svtx_hit_range.first;
0276        hit_iter != svtx_hit_range.second; hit_iter++){
0277
0278     PHG4Hit *this_hit =  hit_iter->second;
0279     assert(this_hit);
0280
0281     e_svtx += this_hit->get_edep();
0282
0283   }
0284
0285   PHG4HitContainer::ConstRange svtx_support_hit_range = _svtx_support_hit_container->getHits();
0286   for (PHG4HitContainer::ConstIterator hit_iter = svtx_support_hit_range.first;
0287        hit_iter != svtx_support_hit_range.second; hit_iter++){
0288
0289     PHG4Hit *this_hit =  hit_iter->second;
0290     assert(this_hit);
0291
0292     e_svtx_support += this_hit->get_edep();
0293
0294   }
0295
0296   //normalize the LCG
0297
0298   lcg = lcg/(e_cemc+ea_cemc+e_emcelect+e_hcin+ea_hcin+e_hcalin_spt+e_hcout+ea_hcout+e_magnet);
0299   hcalin_lcg = hcalin_lcg/(e_hcin+ea_hcin);
0300   hcalout_lcg = hcalout_lcg/(e_hcout+ea_hcout);
0301   emc_lcg = emc_lcg/(e_cemc+ea_cemc);
0302
0303   // Fill the ntuple
0304
0305   float ntdata[27];
0306
0307   ntdata[0] = e_hcin;
0308   ntdata[1] = e_hcout;
0309   ntdata[2] = e_cemc;
0310   ntdata[3] = ea_hcin;
0311   ntdata[4] = ea_hcout;
0312   ntdata[5] = ea_cemc;
0313   ntdata[6] = ev_hcin;
0314   ntdata[7] = ev_hcout;
0315   ntdata[8] = ev_cemc;
0316   ntdata[9] = e_magnet;
0317   ntdata[10] = e_bh;
0318   ntdata[11] = e_emcelect;
0319   ntdata[12] = e_hcalin_spt;
0320
0321   ntdata[13] = e_cemc/(e_cemc + ea_cemc + e_emcelect);
0322   ntdata[14] = e_hcin/(e_hcin + ea_hcin + e_hcalin_spt);
0323   ntdata[15] = e_hcout/(e_hcout + ea_hcout);
0324
0325   ntdata[16] = ev_cemc/(e_cemc + ea_cemc + e_emcelect);
0326   ntdata[17] = ev_hcin/(e_hcin + ea_hcin + e_hcalin_spt);
0327   ntdata[18] = ev_hcout/(e_hcout + ea_hcout);
0328
0329   ntdata[19] = lcg;
0330
0331   ntdata[20] = emc_lcg;
0332   ntdata[21] = hcalin_lcg;
0333   ntdata[22] = hcalout_lcg;
0334
0335   ntdata[23] = e_svtx;
0336   ntdata[24] = e_svtx_support;
0337
0338   ntdata[25] = e_bh_plus;
0339   ntdata[26] = e_bh_minus;
0340
0341   calenergy->Fill(ntdata);
0342
0343   // crosscheck
0344
0345     //std::cout << std::endl;
0346     //std::cout << "ev_cemc = " << ev_cemc << " e_cemc + ea_cemc + e_emcelect = " << e_cemc + ea_cemc + e_emcelect << " ev_cemc/sfvemc " << ntdata[8]/ntdata[16] << std::endl;
0347     //std::cout << "ev_hcin = " << ev_hcin << " e_hcin + ea_hcin + e_hcalin_spt = " << e_hcin + ea_hcin + e_hcalin_spt << " ev_hcin/sfihcal " << ntdata[7]/ntdata[17] << std::endl;
0348     //std::cout << "ev_hcout = " << ev_hcout << " e_hcout + ea_hcout = " << e_hcout + ea_hcout << " ev_hcout/sfohcal " << ntdata[6]/ntdata[18] << std::endl;
0349     //std::cout << std::endl;
0350
0351
0352   // any other return code might lead to aborting the event or analysis for everyone
0353   return 0;
0354 }
0355
0356 int HCALAnalysis::End(PHCompositeNode *topNode)
0357 {
0358
0359   outputfile->Write();
0360   outputfile->Close();
0361
0362   return 0;
0363 }
0364
0365 void HCALAnalysis::GetNodes(PHCompositeNode *topNode)
0366 {
0367
0368   // store top node
0369   _topNode = topNode;
0370
0371   // get DST objects
0372   _hcalout_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALOUT");
0373   if(!_hcalout_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALOUT! No sense continuing" << std::endl; exit(1);}
0374
0375   _hcalin_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALIN");
0376   if(!_hcalin_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALIN! No sense continuing" << std::endl; exit(1);}
0377
0378   _cemc_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_CEMC");
0379   if(!_cemc_hit_container) { std::cout << PHWHERE << ":: No G4HIT_CEMC! No sense continuing" << std::endl; exit(1);}
0380
0381   _hcalout_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_HCALOUT");
0382   if(!_hcalout_abs_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALOUT! No sense continuing" << std::endl; exit(1);}
0383
0384   _hcalin_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_HCALIN");
0385   if(!_hcalin_abs_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALIN! No sense continuing" << std::endl; exit(1);}
0386
0387   _cemc_abs_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_ABSORBER_CEMC");
0388   if(!_cemc_abs_hit_container) { std::cout << PHWHERE << ":: No G4HIT_CEMC! No sense continuing" << std::endl; exit(1);}
0389
0390   _magnet_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_MAGNET");
0391   if(!_magnet_hit_container) { std::cout << PHWHERE << ":: No G4HIT_MAGNET! No sense continuing" << std::endl; exit(1);}
0392
0393   _bh_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_BH_1");
0394   if(!_bh_hit_container) { std::cout << PHWHERE << ":: No G4HIT_BH_1! No sense continuing" << std::endl; exit(1);}
0395
0396   _cemc_electronics_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_CEMC_ELECTRONICS");
0397   if(!_cemc_electronics_hit_container) { std::cout << PHWHERE << ":: No G4HIT_EMCELECTRONICS_0! No sense continuing" << std::endl; exit(1);}
0398
0399   _hcalin_spt_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_HCALIN_SPT");
0400   if(!_hcalin_spt_hit_container) { std::cout << PHWHERE << ":: No G4HIT_HCALIN_SPT! No sense continuing" << std::endl; exit(1);}
0401
0402   _svtx_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_SVTX");
0403   if(!_svtx_hit_container) { std::cout << PHWHERE << ":: No G4HIT_SVTX! No sense continuing" << std::endl; exit(1);}
0404
0405   _svtx_support_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_SVTXSUPPORT");
0406   if(!_svtx_support_hit_container) { std::cout << PHWHERE << ":: No G4HIT_SVTX_SUPPORT! No sense continuing" << std::endl; exit(1);}
0407
0408   _bh_plus_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_BH_FORWARD_PLUS");
0409   if(!_bh_plus_hit_container) { std::cout << PHWHERE << ":: No G4HIT_BH_FORWARD_PLUS! No sense continuing" << std::endl; exit(1);}
0410
0411   _bh_minus_hit_container = findNode::getClass<PHG4HitContainer> (topNode, "G4HIT_BH_FORWARD_NEG");
0412   if(!_bh_minus_hit_container) { std::cout << PHWHERE << ":: No G4HIT_BH_FORWARD_NEG! No sense continuing" << std::endl; exit(1);}
0413
0414   return;
0415 }
0416