Calo_waveform_sim/src/CaloWaveFormSim.cc

0001 #include "CaloWaveFormSim.h"
0002
0003 // G4Hits includes
0004 #include <g4main/PHG4Hit.h>
0005 #include <g4main/PHG4HitContainer.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007 #include <g4main/PHG4Particle.h>
0008 #include <g4main/PHG4HitDefs.h>  // for hit_idbits
0009
0010 // G4Cells includes
0011 #include <g4detectors/PHG4Cell.h>
0012 #include <g4detectors/PHG4CellContainer.h>
0013 #include <g4detectors/PHG4CellDefs.h>  // for genkey, keytype
0014
0015 // Tower includes
0016 #include <calobase/RawTower.h>
0017 #include <calobase/RawTowerContainer.h>
0018 #include <calobase/RawTowerGeom.h>
0019 #include <calobase/RawTowerGeomContainer.h>
0020 #include <calobase/RawTowerDefs.h>
0021 // Cluster includes
0022 #include <calobase/RawCluster.h>
0023 #include <calobase/RawClusterContainer.h>
0024
0025 #include <fun4all/Fun4AllHistoManager.h>
0026 #include <fun4all/Fun4AllReturnCodes.h>
0027
0028
0029 #include "g4detectors/PHG4CylinderGeomContainer.h"
0030 #include "g4detectors/PHG4CylinderGeom_Spacalv1.h"  // for PHG4CylinderGeom_Spaca...
0031 #include "g4detectors/PHG4CylinderGeom_Spacalv3.h"
0032 #include "g4detectors/PHG4CylinderCellGeomContainer.h"
0033 #include "g4detectors/PHG4CylinderCellGeom_Spacalv1.h"
0034
0035
0036 #include <phool/getClass.h>
0037 #include <TProfile.h>
0038 #include <TFile.h>
0039 #include <TNtuple.h>
0040 #include <TMath.h>
0041 #include <cassert>
0042 #include <sstream>
0043 #include <string>
0044 #include <TF1.h>
0045 #include <phool/onnxlib.h>
0046
0047 using namespace std;
0048 TProfile* h_template = nullptr;
0049 TProfile* h_template_ihcal = nullptr;
0050 TProfile* h_template_ohcal = nullptr;
0051
0052
0053 #define ROWDIM 320
0054 #define COLUMNDIM 27
0055
0056
0057
0058 double CaloWaveFormSim::template_function(double *x, double *par)
0059 {
0060   Double_t v1 = par[0]*h_template->Interpolate(x[0]-par[1])+par[2];
0061   return v1;
0062 }
0063
0064 double CaloWaveFormSim::template_function_ihcal(double *x, double *par)
0065 {
0066   Double_t v1 = par[0]*h_template_ihcal->Interpolate(x[0]-par[1])+par[2];
0067   return v1;
0068 }
0069
0070 double CaloWaveFormSim::template_function_ohcal(double *x, double *par)
0071 {
0072   Double_t v1 = par[0]*h_template_ohcal->Interpolate(x[0]-par[1])+par[2];
0073   return v1;
0074 }
0075
0076
0077
0078 CaloWaveFormSim::CaloWaveFormSim(const std::string& name, const std::string& filename)
0079   : SubsysReco(name)
0080   , detector("CEMC")
0081   , outfilename(filename)
0082   , hm(nullptr)
0083   , outfile(nullptr)
0084   , g4hitntuple(nullptr)
0085
0086 {
0087 }
0088
0089 CaloWaveFormSim::~CaloWaveFormSim()
0090 {
0091   delete hm;
0092   delete g4hitntuple;
0093 }
0094
0095 int CaloWaveFormSim::Init(PHCompositeNode*)
0096 {
0097   rnd = new TRandom3(0);
0098   //----------------------------------------------------------------------------------------------------
0099   //Read in the noise file, this currently points to a tim local area file,
0100   //but a copy of this file is in the git repository.
0101   //----------------------------------------------------------------------------------------------------
0102   noise_midrad = new TTree("noise_midrad", "tree");
0103   // noise->ReadFile("/gpfs/mnt/gpfs02/sphenix/user/trinn/fitting_algorithm_playing/slowneutronsignals/CaloAna/noise.csv", "a1:a2:a3:a4:a5:a6:a7:a8:a9:a10:a11:a12:a13:a14:a15:a16:a17:a18:a19:a20:a21:a22:a23:a24:a25:a26:a27:a28:a29:a30:a31");
0104   noise_midrad->ReadFile("/gpfs/mnt/gpfs02/sphenix/user/trinn/sPHENIX_emcal_cosmics_sector0/noise_waveforms/medium_raddmgnoise.csv", "a1:a2:a3:a4:a5:a6:a7:a8:a9:a10:a11:a12:a13:a14:a15:a16:a17:a18:a19:a20:a21:a22:a23:a24:a25:a26:a27:a28:a29:a30:a31");
0105
0106   noise_lowrad = new TTree("noise_lowrad", "tree");
0107   noise_lowrad->ReadFile("/gpfs/mnt/gpfs02/sphenix/user/trinn/sPHENIX_emcal_cosmics_sector0/noise_waveforms/low_raddmgnoise.csv", "a1:a2:a3:a4:a5:a6:a7:a8:a9:a10:a11:a12:a13:a14:a15:a16:a17:a18:a19:a20:a21:a22:a23:a24:a25:a26:a27:a28:a29:a30:a31");
0108
0109
0110   noise_norad = new TTree("noise_norad", "tree");
0111   noise_norad->ReadFile("/gpfs/mnt/gpfs02/sphenix/user/trinn/sPHENIX_emcal_cosmics_sector0/noise_waveforms/no_raddmgnoise.csv", "a1:a2:a3:a4:a5:a6:a7:a8:a9:a10:a11:a12:a13:a14:a15:a16:a17:a18:a19:a20:a21:a22:a23:a24:a25:a26:a27:a28:a29:a30:a31");
0112
0113   for (int i = 0; i < 31;i++)
0114     {
0115       noise_midrad->SetBranchAddress(Form("a%d",i+1),&noise_val_midrad[i]);
0116       noise_lowrad->SetBranchAddress(Form("a%d",i+1),&noise_val_lowrad[i]);
0117       noise_norad->SetBranchAddress(Form("a%d",i+1),&noise_val_norad[i]);
0118     }
0119
0120   hm = new Fun4AllHistoManager(Name());
0121   outfile = new TFile(outfilename.c_str(), "RECREATE");
0122   g4hitntuple = new TTree("tree", "tree");
0123
0124   g4hitntuple->Branch("primpt",&m_primpt);
0125   g4hitntuple->Branch("primeta",&m_primeta);
0126   g4hitntuple->Branch("primphi",&m_primphi);
0127   g4hitntuple->Branch("tedep_emcal",&m_tedep,"tedep_emcal[24576]/F");
0128   g4hitntuple->Branch("tedep_ihcal",&m_tedep_ihcal,"tedep_ihcal[1536]/F");
0129   g4hitntuple->Branch("tedep_ohcal",&m_tedep_ohcal,"tedep_ohcal[1536]/F");
0130   g4hitntuple->Branch("extractedadc_emcal",& m_extractedadc,"extractedadc_emcal[24576]/F");
0131   g4hitntuple->Branch("extractedtime_emcal",& m_extractedtime,"extractedtime_emcal[24576]/F");
0132   g4hitntuple->Branch("extractedadc_ihcal",& m_extractedadc_ihcal,"extractedadc_ihcal[1536]/F");
0133   g4hitntuple->Branch("extractedtime_ihcal",& m_extractedtime_ihcal,"extractedtime_ihcal[1536]/F");
0134   g4hitntuple->Branch("extractedadc_ohcal",& m_extractedadc_ohcal,"extractedadc_ohcal[1536]/F");
0135   g4hitntuple->Branch("extractedtime_ohcal",& m_extractedtime_ohcal,"extractedtime_ohcal[1536]/F");
0136   g4hitntuple->Branch("toweradc_emcal",& m_toweradc,"toweradc_emcal[24576]/F");
0137   g4hitntuple->Branch("toweradc_ihcal",& m_toweradc_ihcal,"toweradc_ihcal[1536]/F");
0138   g4hitntuple->Branch("toweradc_ohcal",& m_toweradc_ohcal,"toweradc_ohcal[1536]/F");
0139   // g4hitntuple->Branch("waveform_ohcal",& m_waveform_ohcal,"waveform_ohcal[1536][16]/I");
0140
0141   // g4hitntuple->Branch("npeaks_ihcal",& m_npeaks_ihcal,"npeaks_ihcall[1536]/I");
0142   // g4hitntuple->Branch("npeaks_ohcal",& m_npeaks_ohcal,"npeaks_ohcal[1536]/I");
0143
0144   //----------------------------------------------------------------------------------------------------
0145   //Read in the template file, this currently points to a tim local area file,
0146   //but a copy of this file is in the git repository.
0147   //----------------------------------------------------------------------------------------------------
0148   // std::string template_input_file = "/gpfs/mnt/gpfs02/sphenix/user/trinn/fitting_algorithm_playing/prdfcode/prototype/offline/packages/Prototype4/templates.root";
0149   std::string cemc_template_input_file = "/gpfs/mnt/gpfs02/sphenix/user/trinn/fitting_algorithm_playing/waveform_simulation/calibrations/WaveformProcessing/templates/testbeam_cemc_template.root";
0150   std::string ihcal_template_input_file = "/gpfs/mnt/gpfs02/sphenix/user/trinn/fitting_algorithm_playing/waveform_simulation/calibrations/WaveformProcessing/templates/testbeam_ihcal_template.root";
0151   std::string ohcal_template_input_file = "/gpfs/mnt/gpfs02/sphenix/user/trinn/fitting_algorithm_playing/waveform_simulation/calibrations/WaveformProcessing/templates/testbeam_ohcal_template.root";
0152
0153   TFile* fin1 = TFile::Open(cemc_template_input_file.c_str());
0154   assert(fin1);
0155   assert(fin1->IsOpen());
0156   h_template = static_cast<TProfile*>(fin1->Get("waveform_template"));
0157
0158   TFile* fin2 = TFile::Open(ihcal_template_input_file.c_str());
0159   assert(fin2);
0160   assert(fin2->IsOpen());
0161
0162   h_template_ihcal = static_cast<TProfile*>(fin2->Get("waveform_template"));
0163
0164   TFile* fin3 = TFile::Open(ohcal_template_input_file.c_str());
0165   assert(fin3);
0166   assert(fin3->IsOpen());
0167
0168   h_template_ohcal = static_cast<TProfile*>(fin3->Get("waveform_template"));
0169
0170
0171   WaveformProcessing = new CaloWaveformProcessing();
0172   WaveformProcessing->set_processing_type(CaloWaveformProcessing::TEMPLATE);
0173   WaveformProcessing->set_template_file("testbeam_cemc_template.root");
0174   WaveformProcessing->initialize_processing();
0175
0176
0177
0178   WaveformProcessing_ihcal = new CaloWaveformProcessing();
0179   WaveformProcessing_ihcal->set_processing_type(CaloWaveformProcessing::TEMPLATE);
0180   WaveformProcessing_ihcal->set_template_file("testbeam_ihcal_template.root");
0181   WaveformProcessing_ihcal->initialize_processing();
0182
0183
0184   WaveformProcessing_ohcal = new CaloWaveformProcessing();
0185   WaveformProcessing_ohcal->set_processing_type(CaloWaveformProcessing::TEMPLATE);
0186   WaveformProcessing_ohcal->set_template_file("testbeam_ohcal_template.root");
0187   WaveformProcessing_ohcal->initialize_processing();
0188
0189
0190   light_collection_model.load_data_file(string(getenv("CALIBRATIONROOT")) + string("/CEMC/LightCollection/Prototype3Module.xml"),
0191                                    "data_grid_light_guide_efficiency", "data_grid_fiber_trans");
0192
0193
0194   std::cout << " hey do you happen? " << std::endl;
0195
0196   return 0;
0197 }
0198
0199 int CaloWaveFormSim::process_event(PHCompositeNode* topNode)
0200 {
0201   process_g4hits(topNode);
0202
0203   return Fun4AllReturnCodes::EVENT_OK;
0204 }
0205
0206 int CaloWaveFormSim::process_g4hits(PHCompositeNode* topNode)
0207 {
0208   bool truth = true;
0209
0210
0211  PHG4CylinderGeomContainer *layergeo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_CEMC");
0212   if (!layergeo)
0213   {
0214     std::cout << "PHG4FullProjSpacalCellReco::process_event - Fatal Error - Could not locate sim geometry node "
0215               <<  std::endl;
0216     exit(1);
0217   }
0218
0219   const PHG4CylinderGeom *layergeom_raw = layergeo->GetFirstLayerGeom();
0220   assert(layergeom_raw);
0221
0222   const PHG4CylinderGeom_Spacalv3 *layergeom =
0223       dynamic_cast<const PHG4CylinderGeom_Spacalv3 *>(layergeom_raw);
0224   assert(layergeom);
0225
0226   PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_CEMC");
0227   PHG4CylinderCellGeom *geo_raw = seggeo->GetFirstLayerCellGeom();
0228   PHG4CylinderCellGeom_Spacalv1 *geo = dynamic_cast<PHG4CylinderCellGeom_Spacalv1 *>(geo_raw);
0229
0230   float waveform[24576][16];
0231   float waveform_ihcal[1536][16];
0232   float waveform_ohcal[1536][16];
0233   float tedep[24576];
0234   float tedep_ihcal[1536];
0235   float tedep_ohcal[1536];
0236   for (int i = 0; i < 24576;i++)
0237     {
0238       for (int j = 0; j < 16;j++)
0239     {
0240       m_waveform[i][j] = 0.0;
0241       waveform[i][j] = 0.0;
0242     }
0243       m_extractedadc[i] = 0;
0244       m_extractedtime[i] = 0;
0245       m_toweradc[i] = 0;
0246       m_ndep[i] = 0.0;
0247       m_tedep[i] = 0.0;
0248       tedep[i] = 0.0;
0249     }
0250
0251
0252   for (int i = 0; i < 1536;i++)
0253     {
0254       tedep_ihcal[i] = 0.0;
0255       m_extractedadc_ihcal[i] = 0;
0256       m_extractedtime_ihcal[i] = 0;
0257       tedep_ohcal[i] = 0.0;
0258       m_extractedadc_ohcal[i] = 0;
0259       m_extractedtime_ohcal[i] = 0;
0260       m_toweradc_ihcal[i] = 0.0;
0261       m_toweradc_ohcal[i] = 0.0;
0262       for (int j = 0; j < 16;j++)
0263     {
0264       m_waveform_ihcal[i][j] = 0;
0265       waveform_ihcal[i][j] = 0;
0266       m_waveform_ohcal[i][j] = 0;
0267       waveform_ohcal[i][j] = 0;
0268     }
0269     }
0270
0271   //---------------------------------------------------------
0272   //Load in the template function as a TF1
0273   //for use in waveform generation
0274   //---------------------------------------------------------
0275
0276   TF1* f_fit = new TF1("f_fit",template_function,0,31,3);
0277   f_fit->SetParameters(1,0,0);
0278
0279   TF1* f_fit_ihcal = new TF1("f_fit_ihcal",template_function_ihcal,0,31,3);
0280   f_fit_ihcal->SetParameters(1,0,0);
0281
0282   TF1* f_fit_ohcal = new TF1("f_fit_ohcal",template_function_ohcal,0,31,3);
0283   f_fit_ohcal->SetParameters(1,0,0);
0284
0285   //-----------------------------------------------------
0286   //Set the timeing in of the prompt
0287   //signal peak to be 4 time samples into
0288   //the waveform
0289   //------------------------------------------------------
0290   float _shiftval = 4-f_fit->GetMaximumX();
0291   f_fit->SetParameters(1,_shiftval,0);
0292
0293   float _shiftval_ihcal = 4-f_fit_ihcal->GetMaximumX();
0294   f_fit_ihcal->SetParameters(1,_shiftval_ihcal,0);
0295
0296   float _shiftval_ohcal = 4-f_fit_ohcal->GetMaximumX();
0297   f_fit_ohcal->SetParameters(1,_shiftval_ohcal,0);
0298
0299
0300   ostringstream nodename;
0301   nodename.str("");
0302   nodename << "G4HIT_" << detector;
0303   PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str().c_str());
0304   PHG4HitContainer* hits_ihcal = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_HCALIN");
0305   PHG4HitContainer* hits_ohcal = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_HCALOUT");
0306
0307   //-----------------------------------------------------------------------------------------
0308   //Loop over truth primary particles and record their kinematics
0309   //-----------------------------------------------------------------------------------------
0310   PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0311   if (truth)
0312     {
0313       PHG4TruthInfoContainer::ConstRange truth_range = truthinfo->GetPrimaryParticleRange();
0314       for (PHG4TruthInfoContainer::ConstIterator truth_iter = truth_range.first; truth_iter != truth_range.second; truth_iter++)
0315     {
0316       PHG4Particle*part = truth_iter->second;
0317       m_primid.push_back(part->get_pid());
0318       float pt =TMath::Sqrt(TMath::Power(part->get_py(),2) +TMath::Power(part->get_px(),2));
0319       float theta = TMath::ATan(pt/part->get_pz());
0320       float phi = TMath::ATan2(part->get_py(),part->get_px());
0321       float eta = -1*TMath::Log(TMath::Tan(fabs(theta/2)));
0322       if (part->get_pz() < 0 )
0323         {
0324           eta = -1 * eta;
0325         }
0326       m_primpt.push_back(part->get_e());
0327       m_primeta.push_back(eta);
0328       m_primphi.push_back(phi);
0329       m_primtrkid.push_back(part->get_track_id());
0330     }
0331     }
0332
0333   if (hits)
0334   {
0335     //-----------------------------------------------------------------------
0336     //Loop over G4Hits to build a waveform simulation
0337     //-----------------------------------------------------------------------
0338     PHG4HitContainer::ConstRange hit_range = hits->getHits();
0339     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0340     {
0341       //-----------------------------------------------------------------
0342       //Extract position information for each G4hit
0343       //information for each G4hit in the calorimeter
0344       //-----------------------------------------------------------------
0345       int scint_id = hit_iter->second->get_scint_id();
0346       PHG4CylinderGeom_Spacalv3::scint_id_coder decoder(scint_id);
0347       std::pair<int, int> tower_z_phi_ID = layergeom->get_tower_z_phi_ID(decoder.tower_ID, decoder.sector_ID);
0348       const int &tower_ID_z = tower_z_phi_ID.first;
0349       const int &tower_ID_phi = tower_z_phi_ID.second;
0350
0351       PHG4CylinderGeom_Spacalv3::tower_map_t::const_iterator it_tower =
0352         layergeom->get_sector_tower_map().find(decoder.tower_ID);
0353       assert(it_tower != layergeom->get_sector_tower_map().end());
0354
0355       const int etabin_cell = geo->get_etabin_block(tower_ID_z);  // block eta bin
0356       const int sub_tower_ID_x = it_tower->second.get_sub_tower_ID_x(decoder.fiber_ID);
0357       const int sub_tower_ID_y = it_tower->second.get_sub_tower_ID_y(decoder.fiber_ID);
0358       unsigned short etabinshort = etabin_cell * layergeom->get_n_subtower_eta() + sub_tower_ID_y;
0359       unsigned short phibin_cell = tower_ID_phi * layergeom->get_n_subtower_phi() + sub_tower_ID_x;
0360
0361       //----------------------------------------------------------------------------------------------------
0362       //Extract light yield from g4hit and correct for light collection efficiency
0363       //----------------------------------------------------------------------------------------------------
0364       double light_yield = hit_iter->second->get_light_yield();
0365       {
0366     const double z = 0.5 * (hit_iter->second->get_local_z(0) + hit_iter->second->get_local_z(1));
0367     assert(not std::isnan(z));
0368     light_yield *= light_collection_model.get_fiber_transmission(z);
0369       }
0370       {
0371     const double x = it_tower->second.get_position_fraction_x_in_sub_tower(decoder.fiber_ID);
0372     const double y = it_tower->second.get_position_fraction_y_in_sub_tower(decoder.fiber_ID);
0373     light_yield *= light_collection_model.get_light_guide_efficiency(x, y);
0374       }
0375       //-------------------------------------------------------------------------
0376       //Map the G4hits to the corresponding CEMC tower
0377       //-------------------------------------------------------------------------
0378       int etabin = etabinshort;
0379       int phibin = phibin_cell;
0380       int towernumber = etabin + 96*phibin;
0381       //---------------------------------------------------------------------
0382       //Convert the G4hit into a waveform contribution
0383       //---------------------------------------------------------------------
0384       // float t0 = 0.5*(hit_iter->second->get_t(0)+hit_iter->second->get_t(1)) / 16.66667;   //Average of g4hit time downscaled by 16.667 ns/time sample
0385       float t0 = (hit_iter->second->get_t(0)) / 16.66667;   //Place waveform at the starting time of the G4hit, avoids issues caused by excessively long lived g4hits
0386       float tmax = 16.667*16;
0387       float tmin = -20;
0388       f_fit->SetParameters(light_yield*26000,_shiftval+t0,0);            //Set the waveform template to match the expected signal from such a hit
0389       if (hit_iter->second->get_t(1) >= tmin && hit_iter->second->get_t(0) <= tmax) {
0390     tedep[towernumber] += light_yield*26000;    // add g4hit adc deposition to the total deposition
0391       }
0392       //-------------------------------------------------------------------------------------------------------------
0393       //For each tower add the new waveform contribution to the total waveform
0394        //-------------------------------------------------------------------------------------------------------------
0395      if (hit_iter->second->get_edep()*26000 > 1 && hit_iter->second->get_t(1) >= tmin && hit_iter->second->get_t(0) <= tmax)
0396     {
0397       for (int j = 0; j < 16;j++) // 16 is the number of time samples
0398         {
0399           waveform[towernumber][j] +=f_fit->Eval(j);
0400         }
0401       m_ndep[towernumber] +=1;
0402     }
0403       m_phibin.push_back(phibin);
0404       m_etabin.push_back(etabin);
0405     }
0406   }
0407
0408
0409
0410   //--------------
0411   // do IHCAL
0412   //--------------
0413
0414   if (hits_ihcal)
0415   {
0416     //-----------------------------------------------------------------------
0417     //Loop over G4Hits to build a waveform simulation
0418     //-----------------------------------------------------------------------
0419     PHG4HitContainer::ConstRange hit_range = hits_ihcal->getHits();
0420     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0421     {
0422       //-----------------------------------------------------------------
0423       //Extract position information for each G4hit
0424       //information for each G4hit in the calorimeter
0425       //-----------------------------------------------------------------
0426       short icolumn = hit_iter->second->get_scint_id();
0427       int introw = (hit_iter->second->get_hit_id() >> PHG4HitDefs::hit_idbits);
0428
0429       if (introw >= ROWDIM || introw < 0)
0430     {
0431       std::cout << __PRETTY_FUNCTION__ << " row " << introw
0432             << " exceed array size: " << ROWDIM
0433             << " adjust ROWDIM and recompile" << std::endl;
0434       exit(1);
0435     }
0436       int towerphi = introw/4;
0437       int towereta = icolumn;
0438
0439       //----------------------------------------------------------------------------------------------------
0440       //Extract light yield from g4hit and correct for light collection efficiency
0441       //----------------------------------------------------------------------------------------------------
0442       double light_yield = hit_iter->second->get_light_yield();  //raw_light_yield has no MEPHI maps applied, light_yield aoppplies the maps change at some point
0443       //-------------------------------------------------------------------------
0444       //Map the G4hits to the corresponding CEMC tower
0445       //-------------------------------------------------------------------------
0446
0447       int etabin = towereta;
0448       int phibin =towerphi;
0449       int towernumber = etabin + 24*phibin;
0450
0451       //---------------------------------------------------------------------
0452       //Convert the G4hit into a waveform contribution
0453       //---------------------------------------------------------------------
0454       // float t0 = 0.5*(hit_iter->second->get_t(0)+hit_iter->second->get_t(1)) / 16.66667;   //Average of g4hit time downscaled by 16.667 ns/time sample
0455       float t0 = (hit_iter->second->get_t(0)) / 16.66667;   //Place waveform at the starting time of the G4hit, avoids issues caused by excessively long lived g4hits
0456       float tmax = 16.667*16;
0457       float tmin = -20;
0458       f_fit_ihcal->SetParameters(light_yield*2600,_shiftval_ihcal+t0,0);            //Set the waveform template to match the expected signal from such a hit
0459       if (hit_iter->second->get_t(1) >= tmin && hit_iter->second->get_t(0) <= tmax) {
0460     tedep_ihcal[towernumber] += light_yield*2600;    // add g4hit adc deposition to the total deposition
0461       }
0462       //-------------------------------------------------------------------------------------------------------------
0463       //For each tower add the new waveform contribution to the total waveform
0464        //-------------------------------------------------------------------------------------------------------------
0465      if (hit_iter->second->get_edep()*2600 > 1 &&hit_iter->second->get_t(1) >= tmin && hit_iter->second->get_t(0) <= tmax )
0466     {
0467       for (int j = 0; j < 16;j++) // 16 is the number of time samples
0468         {
0469           waveform_ihcal[towernumber][j] +=f_fit_ihcal->Eval(j);
0470         }
0471     }
0472     }
0473   }
0474
0475
0476
0477   //--------------
0478   // do OHCAL
0479   //--------------
0480
0481   if (hits_ohcal)
0482   {
0483     //-----------------------------------------------------------------------
0484     //Loop over G4Hits to build a waveform simulation
0485     //-----------------------------------------------------------------------
0486     PHG4HitContainer::ConstRange hit_range = hits_ohcal->getHits();
0487     for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
0488     {
0489       //-----------------------------------------------------------------
0490       //Extract position information for each G4hit
0491       //information for each G4hit in the calorimeter
0492       //-----------------------------------------------------------------
0493       short icolumn = hit_iter->second->get_scint_id();
0494       int introw = (hit_iter->second->get_hit_id() >> PHG4HitDefs::hit_idbits);
0495
0496
0497       if (introw >= ROWDIM || introw < 0)
0498     {
0499       std::cout << __PRETTY_FUNCTION__ << " row " << introw
0500             << " exceed array size: " << ROWDIM
0501             << " adjust ROWDIM and recompile" << std::endl;
0502       exit(1);
0503     }
0504
0505
0506       int towerphi = introw/5;
0507       int towereta = icolumn;
0508
0509       //----------------------------------------------------------------------------------------------------
0510       //Extract light yield from g4hit and correct for light collection efficiency
0511       //----------------------------------------------------------------------------------------------------
0512
0513       double light_yield = hit_iter->second->get_light_yield();  //raw_light_yield has no MEPHI maps applied, light_yield aoppplies the maps change at some point
0514       //-------------------------------------------------------------------------
0515       //Map the G4hits to the corresponding CEMC tower
0516       //-------------------------------------------------------------------------
0517
0518       int etabin = towereta;
0519       int phibin =towerphi;
0520       int towernumber = etabin + 24*phibin;
0521
0522       //---------------------------------------------------------------------
0523       //Convert the G4hit into a waveform contribution
0524       //---------------------------------------------------------------------
0525       // float t0 = 0.5*(hit_iter->second->get_t(0)+hit_iter->second->get_t(1)) / 16.66667;   //Average of g4hit time downscaled by 16.667 ns/time sample
0526       float t0 = (hit_iter->second->get_t(0)) / 16.66667;   //Place waveform at the starting time of the G4hit, avoids issues caused by excessively long lived g4hits
0527       float tmax =16.667*16 ;
0528       float tmin = -20;
0529       f_fit_ohcal->SetParameters(light_yield*5000,_shiftval_ohcal+t0,0);            //Set the waveform template to match the expected signal from such a hit
0530       if (hit_iter->second->get_t(0) < tmax && hit_iter->second->get_t(1) > tmin) {
0531       {
0532         tedep_ohcal[towernumber] += light_yield*5000;    // add g4hit adc deposition to the total deposition
0533       }
0534       }
0535       //-------------------------------------------------------------------------------------------------------------
0536       //For each tower add the new waveform contribution to the total waveform
0537       //-------------------------------------------------------------------------------------------------------------
0538
0539      if (hit_iter->second->get_edep()*5000 > 1 && hit_iter->second->get_t(1) >= tmin && hit_iter->second->get_t(0) <= tmax)
0540     {
0541       for (int j = 0; j < 16;j++) // 16 is the number of time samples
0542         {
0543           waveform_ohcal[towernumber][j] +=f_fit_ohcal->Eval(j);
0544         }
0545     }
0546     }
0547   }
0548
0549
0550   //----------------------------------------------------------------------------------------
0551   //For each tower loop over add a noise waveform
0552   //from cosmic data, gain is too high but is corrected for
0553   //----------------------------------------------------------------------------------------
0554
0555   //-----------------------------
0556   // do noise for EMCal:
0557   //-----------------------------
0558   for (int i = 0; i < 24576;i++)
0559     {
0560       int noise_waveform  = (int)rnd->Uniform(0,1990);
0561       noise_midrad->GetEntry(noise_waveform);
0562       m_tedep[i] = tedep[i];
0563       for (int k = 0; k < 16;k++)
0564     {
0565       m_waveform[i][k] = waveform[i][k]+(noise_val_midrad[k]-1500)/16.0+1500;
0566       // m_waveform[i][k] = waveform[i][k]+1500;
0567     }
0568     }
0569   //---------------------------
0570   // do noise for ihcal:
0571   //---------------------------
0572   for (int i = 0; i < 1536;i++)
0573     {
0574       int noise_waveform  = (int)rnd->Uniform(0,1990);
0575       noise_lowrad->GetEntry(noise_waveform);
0576       m_tedep_ihcal[i] = tedep_ihcal[i];
0577       for (int k = 0; k < 16;k++)
0578     {
0579        m_waveform_ihcal[i][k] = waveform_ihcal[i][k]+(noise_val_lowrad[k]-1500)/16.0+1500;
0580       // m_waveform_ihcal[i][k] = waveform_ihcal[i][k]+1500;
0581     }
0582     }
0583   //---------------------------
0584   // do noise for ohcal:
0585   //---------------------------
0586   for (int i = 0; i < 1536;i++)
0587     {
0588       int noise_waveform  = (int)rnd->Uniform(0,1990);
0589       noise_norad->GetEntry(noise_waveform);
0590
0591       m_tedep_ohcal[i] = tedep_ohcal[i];
0592       for (int k = 0; k < 16;k++)
0593     {
0594       m_waveform_ohcal[i][k] = waveform_ohcal[i][k]+(noise_val_norad[k]-1500)/16.0+1500;
0595       // m_waveform_ohcal[i][k] = waveform_ohcal[i][k]+1500;
0596     }
0597     }
0598   std::vector<std::vector<float>> fitresults;
0599   std::vector<std::vector<float>> fitresults_ihcal;
0600   std::vector<std::vector<float>> fitresults_ohcal;
0601   //------------------------------------------
0602   //Process Waveforms:  EMCal
0603   //------------------------------------------
0604   {
0605     std::vector<std::vector<float>> waveforms;
0606     for (int i = 0; i < 24576;i++)
0607       {
0608     std::vector<float>tmp;
0609     for (int j = 0; j < 16;j++)
0610       {
0611         tmp.push_back(m_waveform[i][j]);
0612       }
0613     waveforms.push_back(tmp);
0614     tmp.clear();
0615       }
0616     fitresults =  WaveformProcessing->process_waveform(waveforms);
0617     for (int i = 0; i < 24576;i++)
0618       {
0619     m_extractedadc[i] = fitresults.at(i).at(0);
0620     m_extractedtime[i] = fitresults.at(i).at(1) - _shiftval;
0621       }
0622     waveforms.clear();
0623   }
0624
0625
0626
0627   //------------------------------------------
0628   //Process Waveforms:  IHCal
0629   //------------------------------------------
0630   {
0631     std::vector<std::vector<float>> waveforms;
0632     for (int i = 0; i < 1536;i++)
0633       {
0634     std::vector<float>tmp;
0635     for (int j = 0; j < 16;j++)
0636       {
0637         tmp.push_back(m_waveform_ihcal[i][j]);
0638       }
0639     waveforms.push_back(tmp);
0640
0641     // int size2 = tmp.size();
0642     // int n_peak = 0;
0643     // for (int j = 2; j < size2-1;j++)
0644     //   {
0645     //     if (tmp.at(j) > 1.01*tmp.at(j-2) && tmp.at(j) > 1.01 * tmp.at(j+1))
0646     //       {
0647     //  n_peak++;
0648     //       }
0649     //   }
0650     // m_npeaks_ihcal[i] = n_peak;
0651
0652
0653
0654     tmp.clear();
0655       }
0656     fitresults_ihcal =  WaveformProcessing_ihcal->process_waveform(waveforms);
0657     for (int i = 0; i < 1536;i++)
0658       {
0659     m_extractedadc_ihcal[i] = fitresults_ihcal.at(i).at(0);
0660     m_extractedtime_ihcal[i] = fitresults_ihcal.at(i).at(1) - _shiftval_ihcal;
0661       }
0662     waveforms.clear();
0663   }
0664
0665
0666
0667   //------------------------------------------
0668   //Process Waveforms:  OHCal
0669   //------------------------------------------
0670   {
0671     std::vector<std::vector<float>> waveforms;
0672     for (int i = 0; i < 1536;i++)
0673       {
0674     std::vector<float>tmp;
0675     for (int j = 0; j < 16;j++)
0676       {
0677         tmp.push_back(m_waveform_ohcal[i][j]);
0678       }
0679     waveforms.push_back(tmp);
0680
0681
0682     // int size2 = tmp.size();
0683     // int n_peak = 0;
0684     // for (int j = 2; j < size2-1;j++)
0685     //   {
0686     //     if (tmp.at(j) - tmp.at(j-2) > 15 && tmp.at(j) - tmp.at(j+1) > 15)
0687     //       {
0688     //  n_peak++;
0689     //       }
0690     //   }
0691
0692     // m_npeaks_ohcal[i] = n_peak;
0693     tmp.clear();
0694       }
0695
0696
0697
0698
0699     fitresults_ohcal =  WaveformProcessing_ohcal->process_waveform(waveforms);
0700     for (int i = 0; i < 1536;i++)
0701       {
0702     m_extractedadc_ohcal[i] = fitresults_ohcal.at(i).at(0);
0703     m_extractedtime_ohcal[i] = fitresults_ohcal.at(i).at(1) - _shiftval_ohcal;
0704       }
0705     waveforms.clear();
0706   }
0707
0708
0709
0710   std::cout << 4 << std::endl;
0711
0712
0713   {
0714     RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_RAW_CEMC");
0715
0716     RawTowerContainer::ConstRange tower_range = towers->getTowers();
0717     for (RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter != tower_range.second; tower_iter++)
0718       {
0719     int phibin = tower_iter->second->get_binphi();
0720     int etabin = tower_iter->second->get_bineta();
0721     float energy = tower_iter->second->get_energy();
0722     int towernumber = etabin + 96*phibin;
0723     m_toweradc[towernumber] =  energy;
0724       }
0725   }
0726
0727   {
0728     RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_RAW_HCALIN");
0729
0730     RawTowerContainer::ConstRange tower_range = towers->getTowers();
0731     for (RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter != tower_range.second; tower_iter++)
0732       {
0733     int phibin = tower_iter->second->get_binphi();
0734     int etabin = tower_iter->second->get_bineta();
0735     float energy = tower_iter->second->get_energy();
0736     int towernumber = etabin + 24*phibin;
0737     m_toweradc_ihcal[towernumber] =  energy;
0738       }
0739   }
0740
0741   {
0742     RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_RAW_HCALOUT");
0743     RawTowerContainer::ConstRange tower_range = towers->getTowers();
0744     for (RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter != tower_range.second; tower_iter++)
0745       {
0746     int phibin = tower_iter->second->get_binphi();
0747     int etabin = tower_iter->second->get_bineta();
0748     float energy = tower_iter->second->get_energy();
0749     int towernumber = etabin + 24*phibin;
0750     m_toweradc_ohcal[towernumber] =  energy;
0751       }
0752   }
0753
0754   g4hitntuple->Fill();
0755
0756
0757   //------------------------------
0758   //Clear vector content
0759   //------------------------------
0760
0761   fitresults.clear();
0762   fitresults_ihcal.clear();
0763   fitresults_ohcal.clear();
0764   m_primid.clear();
0765   m_primtrkid.clear();
0766   m_primpt.clear();
0767   m_primeta.clear();
0768   m_primphi.clear();
0769   m_etabin.clear();
0770   m_phibin.clear();
0771
0772
0773
0774
0775   return Fun4AllReturnCodes::EVENT_OK;
0776 }
0777
0778
0779
0780 int CaloWaveFormSim::End(PHCompositeNode* topNode)
0781 {
0782   outfile->cd();
0783   g4hitntuple->Write();
0784   outfile->Write();
0785   outfile->Close();
0786   delete outfile;
0787   hm->dumpHistos(outfilename, "UPDATE");
0788   return 0;
0789 }