smd_offlmon/src/SmdHistGen.cc

0001 #include "SmdHistGen.h"
0002
0003 #include <fun4all/Fun4AllReturnCodes.h>
0004
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/PHIODataNode.h>    // for PHIODataNode
0007 #include <phool/PHNode.h>          // for PHNode
0008 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0009 #include <phool/PHObject.h>        // for PHObject
0010 #include <phool/getClass.h>
0011
0012 //to analyse PRDF
0013 #include <caloreco/CaloWaveformFitting.h>
0014 #include <ffarawobjects/CaloPacketContainerv1.h>
0015 #include <ffarawobjects/CaloPacket.h>
0016 /* #include <Event/Event.h> */
0017 /* #include <Event/EventTypes.h> */
0018 /* #include <Event/packet.h> */
0019
0020 //to analyse DST
0021 #include <calobase/TowerInfoContainer.h>
0022 #include <calobase/TowerInfo.h>
0023 #include <ffaobjects/RunHeaderv1.h>
0024 #include <ffarawobjects/Gl1Packetv1.h>
0025
0026 //spin database stuff
0027 #include <uspin/SpinDBContent.h>
0028 #include <uspin/SpinDBOutput.h>
0029
0030 //ROOT stuff
0031 #include <TFile.h>
0032 #include <TH1.h>
0033 #include <TH2.h>
0034 /* #include <TGraphErrors.h> */
0035 #include <TString.h>
0036 #include <TCanvas.h>
0037 #include <TLegend.h>
0038 #include <TStyle.h>
0039
0040 #include <algorithm>
0041
0042 //____________________________________________________________________________..
0043 SmdHistGen::SmdHistGen(const std::string& name, const std::string& which_mode, const char* outname):
0044   SubsysReco(name),
0045   mode(which_mode),
0046   outfilename(outname),
0047   outfile(nullptr)
0048 {
0049   std::cout << "SmdHistGen::SmdHistGen(const std::string &name) Calling ctor" << std::endl;
0050 }
0051
0052 //____________________________________________________________________________..
0053 SmdHistGen::~SmdHistGen()
0054 {
0055   std::cout << "SmdHistGen::~SmdHistGen() Calling dtor" << std::endl;
0056 }
0057
0058 //____________________________________________________________________________..
0059 int SmdHistGen::Init(PHCompositeNode *topNode)
0060 {
0061   std::cout << "SmdHistGen::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0062
0063   // Create output file
0064   outfile = new TFile(outfilename, "RECREATE");
0065   outfile->cd();
0066
0067   // Create hitograms
0068
0069   // hits per channel and hit multiply
0070   smd_hor_north_total_hits = new TH1F("smd_hor_north_total_hits", "Total Hits by Channel, SMD North Horizontal;Channel;Counts", 8, -0.5, 7.5);
0071   smd_ver_north_total_hits = new TH1F("smd_ver_north_total_hits", "Total Hits by Channel, SMD North Vertical;Channel;Counts", 7, -0.5, 6.5);
0072   smd_hor_north_neutron_hits = new TH1F("smd_hor_north_neutron_hits", "Neutron Hits by Channel, SMD North Horizontal;Channel;Counts", 8, -0.5, 7.5);
0073   smd_ver_north_neutron_hits = new TH1F("smd_ver_north_neutron_hits", "Neutron Hits by Channel, SMD North Vertical;Channel;Counts", 7, -0.5, 6.5);
0074   smd_hor_north_multiplicity = new TH1F("smd_hor_north_multiplicity", Form("In-Time Hit Multiplicity, ADC>%d, SMD North Horizontal;# Hits;Counts", minSMDcut), 9, -0.5, 8.5);
0075   smd_ver_north_multiplicity = new TH1F("smd_ver_north_multiplicity", Form("In-Time Hit Multiplicity, ADC>%d, SMD North Vertical;# Hits;Counts", minSMDcut), 8, -0.5, 7.5);
0076   smd_hor_south_total_hits = new TH1F("smd_hor_south_total_hits", "Total Hits by Channel, SMD South Horizontal;Channel;Counts", 8, -0.5, 7.5);
0077   smd_ver_south_total_hits = new TH1F("smd_ver_south_total_hits", "Total Hits by Channel, SMD South Vertical;Channel;Counts", 7, -0.5, 6.5);
0078   smd_hor_south_neutron_hits = new TH1F("smd_hor_south_neutron_hits", "Neutron Hits by Channel, SMD South Horizontal;Channel;Counts", 8, -0.5, 7.5);
0079   smd_ver_south_neutron_hits = new TH1F("smd_ver_south_neutron_hits", "Neutron Hits by Channel, SMD South Vertical;Channel;Counts", 7, -0.5, 6.5);
0080   smd_hor_south_multiplicity = new TH1F("smd_hor_south_multiplicity", Form("In-Time Hit Multiplicity, ADC>%d, SMD South Horizontal;# Hits;Counts", minSMDcut), 9, -0.5, 8.5);
0081   smd_ver_south_multiplicity = new TH1F("smd_ver_south_multiplicity", Form("In-Time Hit Multiplicity, ADC>%d, SMD South Vertical;# Hits;Counts", minSMDcut), 8, -0.5, 7.5);
0082
0083   // signals and waveforms
0084   float signal_max = 5000;
0085   smd_north_signal = new TH1F("smd_north_signal", "North SMD Peak ADC, All Channels;ADC;Counts", 1000, 0.0, signal_max);
0086   smd_south_signal = new TH1F("smd_south_signal", "South SMD Peak ADC, All Channels;ADC;Counts", 1000, 0.0, signal_max);
0087   for (int i=0; i<15; i++)
0088   {
0089     smd_north_signals[i] = new TH1F(Form("smd_north_signal_%d", i), Form("North SMD Channel %d", i), 1000, 0.0, signal_max);
0090     smd_north_channel_waveforms[i] = new TH2F(Form("smd_north_waveforms_%d", i), Form("North SMD Waveforms Channel %d", i), 128, -0.5, 16.5, 256, 0.0, signal_max);
0091     smd_south_signals[i] = new TH1F(Form("smd_south_signal_%d", i), Form("South SMD Channel %d", i), 1000, 0.0, signal_max);
0092     smd_south_channel_waveforms[i] = new TH2F(Form("smd_south_waveforms_%d", i), Form("South SMD Waveforms Channel %d", i), 128, -0.5, 16.5, 256, 0.0, signal_max);
0093   }
0094   smd_sum_hor_north = new TH1F ("smd_sum_hor_north", "SMD North y", 512, 0, 2048);
0095   smd_sum_ver_north = new TH1F ("smd_sum_ver_north", "SMD North x", 512, 0, 2048);
0096   smd_sum_hor_south = new TH1F ("smd_sum_hor_south", "SMD South y", 640, 0, 2560);
0097   smd_sum_ver_south = new TH1F ("smd_sum_ver_south", "SMD South x", 640, 0, 2560);
0098   smd_north_waveforms = new TH2F("smd_north_waveforms", "North SMD Waveform;Time;Max ADC", 128, -0.5, 16.5, 256, 0.0, signal_max);
0099   smd_north_waveforms_zoomed = new TH2F("smd_north_waveforms_zoomed", "North SMD Waveform;Time;Max ADC", 128, -0.5, 16.5, 250, -0.5, 249.5);
0100   smd_north_waveforms_raw = new TH2F("smd_north_waveforms_raw", "North SMD Waveform;Time;Max ADC", 17, -0.5, 16.5, 1000, 0.0, signal_max);
0101   smd_south_waveforms = new TH2F("smd_south_waveforms", "South SMD Waveform;Time;Max ADC", 128, -0.5, 16.5, 256, 0.0, signal_max);
0102   smd_south_waveforms_zoomed = new TH2F("smd_south_waveforms_zoomed", "South SMD Waveform;Time;Max ADC", 128, -0.5, 16.5, 250, -0.5, 249.5);
0103   smd_south_waveforms_raw = new TH2F("smd_south_waveforms_raw", "South SMD Waveform;Time;Max ADC", 17, -0.5, 16.5, 1000, 0.0, signal_max);
0104   smd_north_pedestal = new TH1F("smd_north_pedestal", "North SMD Pedestal;ADC;Counts", 1000, 0.0, signal_max);
0105   smd_south_pedestal = new TH1F("smd_south_pedestal", "South SMD Pedestal;ADC;Counts", 1000, 0.0, signal_max);
0106   zdc1_north = new TH1F("zdc1_north", "North ZDC1 Signal;ADC;Counts", 256, 0, signal_max);
0107   zdc2_north = new TH1F("zdc2_north", "North ZDC2 Signal;ADC;Counts", 256, 0, signal_max);
0108   zdc_north_waveforms = new TH2F("zdc_north_waveforms", "North ZDC Waveform;Time;Max ADC", 128, -0.5, 16.5, 256, 0.0, signal_max);
0109   zdc1_south = new TH1F("zdc1_south", "South ZDC1 Signal;ADC;Counts", 256, 0, signal_max);
0110   zdc2_south = new TH1F("zdc2_south", "South ZDC2 Signal;ADC;Counts", 256, 0, signal_max);
0111   zdc_south_waveforms = new TH2F("zdc_south_waveforms", "South ZDC Waveform;Time;Max ADC", 128, -0.5, 16.5, 256, 0.0, signal_max);
0112   vetofront_north = new TH1F("vetofront_north", "North Front Veto Signal;ADC;Counts", 256, 0, signal_max);
0113   vetoback_north = new TH1F("vetoback_north", "North Back Veto Signal;ADC;Counts", 256, 0, signal_max);
0114   veto_north_waveforms = new TH2F("veto_north_waveforms", "North Veto Waveform;Time;Max ADC", 128, -0.5, 16.5, 256, 0.0, signal_max);
0115   vetofront_south = new TH1F("vetofront_south", "South Front Veto Signal;ADC;Counts", 256, 0, signal_max);
0116   vetoback_south = new TH1F("vetoback_south", "South Back Veto Signal;ADC;Counts", 256, 0, signal_max);
0117   veto_south_waveforms = new TH2F("veto_south_waveforms", "South Veto Waveform;Time;Max ADC", 128, -0.5, 16.5, 256, 0.0, signal_max);
0118
0119   // xy distributions
0120   int nbins_xy = 50;
0121   smd_xy_all_north = new TH2F("smd_xy_all_north", "SMD hit position north, all good hits;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0122   smd_xy_all_corrected_north = new TH2F("smd_xy_all_corrected_north", "Center-Corrected SMD hit position north, all good hits;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0123   smd_xy_neutron_north = new TH2F("smd_xy_neutron_north", "SMD hit position north, neutron hits only;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0124   smd_xy_neutron_corrected_north = new TH2F("smd_xy_neutron_corrected_north", "Center-Corrected SMD hit position north, neutron hits only;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0125   smd_xy_all_south = new TH2F("smd_xy_all_south", "SMD hit position south, all good hits;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0126   smd_xy_all_corrected_south = new TH2F("smd_xy_all_corrected_south", "Center-Corrected SMD hit position south, all good hits;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0127   smd_xy_neutron_south = new TH2F("smd_xy_neutron_south", "SMD hit position south, neutron hits only;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0128   smd_xy_neutron_corrected_south = new TH2F("smd_xy_neutron_corrected_south", "Center-Corrected SMD hit position south, neutron hits only;x;y", 110, -5.5, 5.5, 119, -5.92, 5.92);
0129
0130   // 1D spin dependent x and y
0131   smd_hor_north = new TH1F("smd_hor_north", "Neutron centroid distribution, SMD North y", nbins_xy, -5.92, 5.92);
0132   smd_ver_north = new TH1F("smd_ver_north", "Neutron centroid distribution, SMD North x", nbins_xy, -5.5, 5.5);
0133   smd_hor_north_up = new TH1F("smd_hor_north_up", "Neutron centroid distribution, SMD North y, Spin Up", nbins_xy, -5.92, 5.92);
0134   smd_ver_north_up = new TH1F("smd_ver_north_up", "Neutron centroid distribution, SMD North x, Spin Up", nbins_xy, -5.5, 5.5);
0135   smd_hor_north_down = new TH1F("smd_hor_north_down", "Neutron centroid distribution, SMD North y, Spin Down", nbins_xy, -5.92, 5.92);
0136   smd_ver_north_down = new TH1F("smd_ver_north_down", "Neutron centroid distribution, SMD North x, Spin Down", nbins_xy, -5.5, 5.5);
0137   smd_hor_south = new TH1F("smd_hor_south", "Neutron centroid distribution, SMD South y", nbins_xy, -5.92, 5.92);
0138   smd_ver_south = new TH1F("smd_ver_south", "Neutron centroid distribution, SMD South x", nbins_xy, -5.5, 5.5);
0139   smd_hor_south_up = new TH1F("smd_hor_south_up", "Neutron centroid distribution, SMD South y, Spin Up", nbins_xy, -5.92, 5.92);
0140   smd_ver_south_up = new TH1F("smd_ver_south_up", "Neutron centroid distribution, SMD South x, Spin Up", nbins_xy, -5.5, 5.5);
0141   smd_hor_south_down = new TH1F("smd_hor_south_down", "Neutron centroid distribution, SMD South y, Spin Down", nbins_xy, -5.92, 5.92);
0142   smd_ver_south_down = new TH1F("smd_ver_south_down", "Neutron centroid distribution, SMD South x, Spin Down", nbins_xy, -5.5, 5.5);
0143
0144   // beam center correction
0145   smd_hor_north_corrected = new TH1F("smd_hor_north_corrected", "Center-Corrected Neutron centroid distribution, SMD North y", nbins_xy, -5.92, 5.92);
0146   smd_ver_north_corrected = new TH1F("smd_ver_north_corrected", "Center-Corrected Neutron centroid distribution, SMD North x", nbins_xy, -5.5, 5.5);
0147   smd_hor_north_corrected_up = new TH1F("smd_hor_north_corrected_up", "Center-Corrected Neutron centroid distribution, SMD North y, Spin Up", nbins_xy, -5.92, 5.92);
0148   smd_ver_north_corrected_up = new TH1F("smd_ver_north_corrected_up", "Center-Corrected Neutron centroid distribution, SMD North x, Spin Up", nbins_xy, -5.5, 5.5);
0149   smd_hor_north_corrected_down = new TH1F("smd_hor_north_corrected_down", "Center-Corrected Neutron centroid distribution, SMD North y, Spin Down", nbins_xy, -5.92, 5.92);
0150   smd_ver_north_corrected_down = new TH1F("smd_ver_north_corrected_down", "Center-Corrected Neutron centroid distribution, SMD North x, Spin Down", nbins_xy, -5.5, 5.5);
0151   smd_hor_south_corrected = new TH1F("smd_hor_south_corrected", "Center-Corrected Neutron centroid distribution, SMD South y", nbins_xy, -5.92, 5.92);
0152   smd_ver_south_corrected = new TH1F("smd_ver_south_corrected", "Center-Corrected Neutron centroid distribution, SMD South x", nbins_xy, -5.5, 5.5);
0153   smd_hor_south_corrected_up = new TH1F("smd_hor_south_corrected_up", "Center-Corrected Neutron centroid distribution, SMD South y, Spin Up", nbins_xy, -5.92, 5.92);
0154   smd_ver_south_corrected_up = new TH1F("smd_ver_south_corrected_up", "Center-Corrected Neutron centroid distribution, SMD South x, Spin Up", nbins_xy, -5.5, 5.5);
0155   smd_hor_south_corrected_down = new TH1F("smd_hor_south_corrected_down", "Center-Corrected Neutron centroid distribution, SMD South y, Spin Down", nbins_xy, -5.92, 5.92);
0156   smd_ver_south_corrected_down = new TH1F("smd_ver_south_corrected_down", "Center-Corrected Neutron centroid distribution, SMD South x, Spin Down", nbins_xy, -5.5, 5.5);
0157
0158   // phi distributions for asymmetry
0159   smd_north_phi_up = new TH1F("smd_north_phi_up", "Spin-up #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0160   smd_north_phi_down = new TH1F("smd_north_phi_down", "Spin-down #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0161   smd_north_phi_sum = new TH1F("smd_north_phi_sum", "Up+Down #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0162   smd_north_phi_diff = new TH1F("smd_north_phi_diff", "Up-Down #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0163   smd_north_phi_L_up = new TH1F("smd_north_phi_L_up", "Spin-up #phi distribution, SMD North;#phi;Counts", 8, -PI, 0.0);
0164   smd_north_phi_L_down = new TH1F("smd_north_phi_L_down", "Spin-down #phi distribution, SMD North;#phi;Counts", 8, -PI, 0.0);
0165   smd_north_phi_R_up = new TH1F("smd_north_phi_R_up", "Spin-up #phi distribution, SMD North;#phi;Counts", 8, 0.0, PI);
0166   smd_north_phi_R_down = new TH1F("smd_north_phi_R_down", "Spin-down #phi distribution, SMD North;#phi;Counts", 8, 0.0, PI);
0167   smd_south_phi_up = new TH1F("smd_south_phi_up", "Spin-up #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0168   smd_south_phi_down = new TH1F("smd_south_phi_down", "Spin-down #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0169   smd_south_phi_sum = new TH1F("smd_south_phi_sum", "Up+Down #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0170   smd_south_phi_diff = new TH1F("smd_south_phi_diff", "Up-Down #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0171   smd_south_phi_L_up = new TH1F("smd_south_phi_L_up", "Spin-up #phi distribution, SMD South;#phi;Counts", 8, -PI, 0.0);
0172   smd_south_phi_L_down = new TH1F("smd_south_phi_L_down", "Spin-down #phi distribution, SMD South;#phi;Counts", 8, -PI, 0.0);
0173   smd_south_phi_R_up = new TH1F("smd_south_phi_R_up", "Spin-up #phi distribution, SMD South;#phi;Counts", 8, 0.0, PI);
0174   smd_south_phi_R_down = new TH1F("smd_south_phi_R_down", "Spin-down #phi distribution, SMD South;#phi;Counts", 8, 0.0, PI);
0175
0176   // beam center correction
0177   smd_north_phi_up_corrected = new TH1F("smd_north_phi_up_corrected", "Center-Corrected Spin-up #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0178   smd_north_phi_down_corrected = new TH1F("smd_north_phi_down_corrected", "Center-Corrected Spin-down #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0179   smd_north_phi_sum_corrected = new TH1F("smd_north_phi_sum_corrected", "Center-Corrected Up+Down #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0180   smd_north_phi_diff_corrected = new TH1F("smd_north_phi_diff_corrected", "Center-Corrected Up-Down #phi distribution, SMD North;#phi;Counts", 8, -PI, PI);
0181   smd_north_phi_L_up_corrected = new TH1F("smd_north_phi_L_up_corrected", "Center-Corrected Spin-up #phi distribution, SMD North;#phi;Counts", 8, -PI, 0.0);
0182   smd_north_phi_L_down_corrected = new TH1F("smd_north_phi_L_down_corrected", "Center-Corrected Spin-down #phi distribution, SMD North;#phi;Counts", 8, -PI, 0.0);
0183   smd_north_phi_R_up_corrected = new TH1F("smd_north_phi_R_up_corrected", "Center-Corrected Spin-up #phi distribution, SMD North;#phi;Counts", 8, 0.0, PI);
0184   smd_north_phi_R_down_corrected = new TH1F("smd_north_phi_R_down_corrected", "Center-Corrected Spin-down #phi distribution, SMD North;#phi;Counts", 8, 0.0, PI);
0185   smd_south_phi_up_corrected = new TH1F("smd_south_phi_up_corrected", "Center-Corrected Spin-up #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0186   smd_south_phi_down_corrected = new TH1F("smd_south_phi_down_corrected", "Center-Corrected Spin-down #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0187   smd_south_phi_sum_corrected = new TH1F("smd_south_phi_sum_corrected", "Center-Corrected Up+Down #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0188   smd_south_phi_diff_corrected = new TH1F("smd_south_phi_diff_corrected", "Center-Corrected Up-Down #phi distribution, SMD South;#phi;Counts", 8, -PI, PI);
0189   smd_south_phi_L_up_corrected = new TH1F("smd_south_phi_L_up_corrected", "Center-Corrected Spin-up #phi distribution, SMD South;#phi;Counts", 8, -PI, 0.0);
0190   smd_south_phi_L_down_corrected = new TH1F("smd_south_phi_L_down_corrected", "Center-Corrected Spin-down #phi distribution, SMD South;#phi;Counts", 8, -PI, 0.0);
0191   smd_south_phi_R_up_corrected = new TH1F("smd_south_phi_R_up_corrected", "Center-Corrected Spin-up #phi distribution, SMD South;#phi;Counts", 8, 0.0, PI);
0192   smd_south_phi_R_down_corrected = new TH1F("smd_south_phi_R_down_corrected", "Center-Corrected Spin-down #phi distribution, SMD South;#phi;Counts", 8, 0.0, PI);
0193
0194   // Set relative gain values to unity for now
0195   for (int i=0; i<16; i++) {
0196     smd_north_rgain[i] = 1.0;
0197     smd_south_rgain[i] = 1.0;
0198   }
0199
0200   // The next block will be how we read the relative gain values from a
0201   // calibrations file
0202   /*
0203   // read our calibrations from ZdcMonData.dat
0204   const char *zdccalib = getenv("ZDCCALIB");
0205   if (!zdccalib)
0206   {
0207     std::cout << "ZDCCALIB environment variable not set" << std::endl;
0208     exit(1);
0209   }
0210
0211   //getting gains
0212   float col1, col2, col3;
0213   std::string gainfile = std::string(zdccalib) + "/" + "/ZdcCalib.pmtgain";
0214   std::ifstream gain_infile(gainfile);
0215
0216   if (!gain_infile)
0217   {
0218     std::cout << gainfile << " could not be opened." ;
0219     exit(1);
0220   }
0221
0222   for (int i = 0; i < 32; i++)
0223   {
0224     gain_infile >> col1 >> col2 >> col3;
0225     gain[i] = col1;
0226   }
0227
0228   for (int i = 0; i < 16; i++)  // relative gains of SMD channels
0229   {
0230     smd_south_rgain[i] = gain[i];  // 0-7: y channels, 8-14: x channels, 15: analog sum
0231     smd_north_rgain[i] = gain[i + 16];  // same as above
0232   }
0233
0234   gain_infile.close();
0235   */
0236   // done reading gains from file
0237
0238   return Fun4AllReturnCodes::EVENT_OK;
0239 }
0240
0241 //____________________________________________________________________________..
0242 int SmdHistGen::InitRun(PHCompositeNode *topNode)
0243 {
0244   std::cout << "SmdHistGen::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0245
0246   // Get run number
0247   runHeader = findNode::getClass<RunHeaderv1>(topNode, "RunHeader");
0248   if (!runHeader)
0249   {
0250     std::cout << PHWHERE << ":: Missing RunHeader! Exiting!" << std::endl;
0251     exit(1);
0252   }
0253   runNum = runHeader->get_RunNumber();
0254   std::cout << "Run number is " << runNum << std::endl;
0255
0256   // Get spin pattern info
0257   int spinDB_status = GetSpinPatterns();
0258   if (spinDB_status)
0259   {
0260     std::cout << PHWHERE << ":: Run number " << runNum << " not found in spin DB! Exiting!" << std::endl;
0261     /* exit(1); */
0262   }
0263
0264   return Fun4AllReturnCodes::EVENT_OK;
0265 }
0266
0267 //____________________________________________________________________________..
0268 int SmdHistGen::process_event(PHCompositeNode *topNode)
0269 {
0270   /* std::cout << "SmdHistGen::process_event(PHCompositeNode *topNode) Processing Event" << std::endl; */
0271
0272   evtctr++;
0273   if (evtctr%10000 == 0) std::cout << "Event " << evtctr << std::endl;
0274
0275   // Get the bunch number for this event
0276   Gl1Packetv1* gl1 = findNode::getClass<Gl1Packetv1>(topNode, "GL1Packet");
0277   if (gl1)
0278   {
0279     bunchNum = gl1->getBunchNumber();
0280     bunchNum = (bunchNum + crossingShift)%NBUNCHES;
0281     /* std::cout << "Got bunch number = " << bunchNum << std::endl; */
0282   }
0283   else
0284   {
0285     // for testing
0286     std::cout << "GL1 missing!" << std::endl;
0287     bunchNum = evtctr%4;
0288   }
0289   /* gl1event++; */
0290   /* std::cout << Form("evtctr=%d, gl1event=%d, smdevent=%d", evtctr, gl1event, smdevent) << std::endl; */
0291
0292   // To account for the GL1-ZDC event # mismatch:
0293   // Compute the position and asymmetry-related parts *before* getting the new
0294   // ADC values for this event -- that way, the hit positions etc. are computed
0295   // based on the *previous* event's ADC values
0296
0297   /* if (evtctr == 1) */
0298   /* { */
0299   /*   towerinfosZDC = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_ZDC"); */
0300   /*   packetsZDC = findNode::getClass<CaloPacketContainerv1>(topNode, "ZDCPackets"); */
0301   /*   GetAdcs(); */
0302   /*   return Fun4AllReturnCodes::EVENT_OK; */
0303   /* } */
0304
0305   /*
0306   // Get the ADC values
0307   towerinfosZDC = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_ZDC");
0308   packetsZDC = findNode::getClass<CaloPacketContainerv1>(topNode, "ZDCPackets");
0309   GetAdcs();
0310   */
0311
0312   // call the functions
0313   CompSmdAdc();
0314   CompSmdPos();
0315   /* std::cout << Form("Event %d start: smd_adc = {", evtctr); */
0316   /* for (int i=0; i<15; i++) { */
0317   /*     std::cout << smd_adc[i] << ", "; */
0318   /* } */
0319   /* std::cout << smd_adc[15] << "}" << std::endl; */
0320   /* std::cout << Form("Event %d start: smd_time = {", evtctr); */
0321   /* for (int i=0; i<15; i++) { */
0322   /*     std::cout << smd_time[i] << ", "; */
0323   /* } */
0324   /* std::cout << smd_time[15] << "}" << std::endl; */
0325   CompSmdPosCorr();
0326   CompSumSmd();
0327   CountSMDHits();
0328   if ((n_hor_numhits > minSMDhits) && (n_ver_numhits > minSMDhits))
0329   {
0330     smd_xy_all_north->Fill(smd_pos[1], smd_pos[0]);
0331     smd_xy_all_corrected_north->Fill(smd_pos_corr[1], smd_pos_corr[0]);
0332   }
0333   if ((s_hor_numhits > minSMDhits) && (s_ver_numhits > minSMDhits))
0334   {
0335     smd_xy_all_south->Fill(smd_pos[3], smd_pos[2]);
0336     smd_xy_all_corrected_south->Fill(smd_pos_corr[3], smd_pos_corr[2]);
0337   }
0338   bool n_neutron = NeutronSelection("north");
0339   bool s_neutron = NeutronSelection("south");
0340
0341   // FILLING OUT THE HISTOGRAMS
0342   if (n_neutron)
0343   {
0344     smd_hor_north->Fill(smd_pos[0]);
0345     smd_ver_north->Fill(smd_pos[1]);
0346     smd_sum_ver_north->Fill(smd_sum[1]);
0347     smd_sum_hor_north->Fill(smd_sum[0]);
0348     smd_xy_neutron_north->Fill(smd_pos[1], smd_pos[0]);
0349     //std::cout<<" smd sum ver north = "<<smd_sum[1]<<std::endl;
0350     //std::cout<<" smd sum hor north = "<<smd_sum[0]<<std::endl;
0351
0352     // SMD hits by channel
0353     for ( int i = 0; i < 8; i++)
0354     {
0355       if ( smd_adc[i] > minSMDcut ) {smd_hor_north_neutron_hits->Fill(i);}
0356     }
0357     for ( int i = 0; i < 7; i++)
0358     {
0359       if ( smd_adc[i + 8] > minSMDcut ) {smd_ver_north_neutron_hits->Fill(i);}
0360     }
0361
0362     // Separate spin up and down centroid positions
0363     float north_x = smd_pos[1];
0364     float north_y = smd_pos[0];
0365     float north_phi = atan2(north_y, north_x);
0366     int blueSpin = spinPatternBlue[bunchNum];
0367     if (blueSpin == 1)
0368     {
0369       smd_hor_north_up->Fill(north_y);
0370       smd_ver_north_up->Fill(north_x);
0371       smd_north_phi_up->Fill(north_phi);
0372       smd_north_phi_L_up->Fill(north_phi);
0373       smd_north_phi_R_up->Fill(north_phi);
0374     }
0375     if (blueSpin == -1)
0376     {
0377       smd_hor_north_down->Fill(north_y);
0378       smd_ver_north_down->Fill(north_x);
0379       smd_north_phi_down->Fill(north_phi);
0380       smd_north_phi_L_down->Fill(north_phi);
0381       smd_north_phi_R_down->Fill(north_phi);
0382     }
0383   }
0384
0385   bool n_neutron_corr = NeutronSelection("north", true);
0386   if (n_neutron_corr)
0387   {
0388     smd_hor_north_corrected->Fill(smd_pos_corr[0]);
0389     smd_ver_north_corrected->Fill(smd_pos_corr[1]);
0390     smd_xy_neutron_corrected_north->Fill(smd_pos_corr[1], smd_pos_corr[0]);
0391
0392     // Separate spin up and down centroid positions
0393     float north_x_corr = smd_pos_corr[1];
0394     float north_y_corr = smd_pos_corr[0];
0395     float north_phi_corr = atan2(north_y_corr, north_x_corr);
0396     int blueSpin = spinPatternBlue[bunchNum];
0397     if (blueSpin == 1)
0398     {
0399       smd_hor_north_corrected_up->Fill(north_y_corr);
0400       smd_ver_north_corrected_up->Fill(north_x_corr);
0401       smd_north_phi_up_corrected->Fill(north_phi_corr);
0402       smd_north_phi_L_up_corrected->Fill(north_phi_corr);
0403       smd_north_phi_R_up_corrected->Fill(north_phi_corr);
0404     }
0405     if (blueSpin == -1)
0406     {
0407       smd_hor_north_corrected_down->Fill(north_y_corr);
0408       smd_ver_north_corrected_down->Fill(north_x_corr);
0409       smd_north_phi_down_corrected->Fill(north_phi_corr);
0410       smd_north_phi_L_down_corrected->Fill(north_phi_corr);
0411       smd_north_phi_R_down_corrected->Fill(north_phi_corr);
0412     }
0413   }
0414
0415   if (s_neutron)
0416   {
0417     smd_hor_south->Fill(smd_pos[2]);
0418     smd_ver_south->Fill(smd_pos[3]);
0419     smd_sum_ver_south->Fill(smd_sum[3]);
0420     smd_sum_hor_south->Fill(smd_sum[2]);
0421     smd_xy_neutron_south->Fill(smd_pos[3], smd_pos[2]);
0422     //std::cout<<" smd sum hor south = "<<smd_sum[2]<<std::endl;
0423     //std::cout<<" smd sum ver south = "<<smd_sum[3]<<std::endl;
0424
0425     // SMD hits by channel
0426     for ( int i = 0; i < 8; i++)
0427     {
0428       if ( smd_adc[i] > minSMDcut ) {smd_hor_south_neutron_hits->Fill(i);}
0429     }
0430     for ( int i = 0; i < 7; i++)
0431     {
0432       if ( smd_adc[i + 8] > minSMDcut ) {smd_ver_south_neutron_hits->Fill(i);}
0433     }
0434
0435     // Separate spin up and down centroid positions
0436     float south_x = smd_pos[3];
0437     float south_y = smd_pos[2];
0438     float south_phi = atan2(south_y, south_x);
0439     int yellowSpin = spinPatternYellow[bunchNum];
0440     if (yellowSpin == 1)
0441     {
0442       smd_hor_south_up->Fill(south_y);
0443       smd_ver_south_up->Fill(south_x);
0444       smd_south_phi_up->Fill(south_phi);
0445       smd_south_phi_L_up->Fill(south_phi);
0446       smd_south_phi_R_up->Fill(south_phi);
0447     }
0448     if (yellowSpin == -1)
0449     {
0450       smd_hor_south_down->Fill(south_y);
0451       smd_ver_south_down->Fill(south_x);
0452       smd_south_phi_down->Fill(south_phi);
0453       smd_south_phi_L_down->Fill(south_phi);
0454       smd_south_phi_R_down->Fill(south_phi);
0455     }
0456   }
0457
0458   bool s_neutron_corr = NeutronSelection("south", true);
0459   if (s_neutron_corr)
0460   {
0461     smd_hor_south_corrected->Fill(smd_pos_corr[2]);
0462     smd_ver_south_corrected->Fill(smd_pos_corr[3]);
0463     smd_xy_neutron_corrected_south->Fill(smd_pos_corr[3], smd_pos_corr[2]);
0464     //std::cout<<" smd sum hor south = "<<smd_sum[2]<<std::endl;
0465     //std::cout<<" smd sum ver south = "<<smd_sum[3]<<std::endl;
0466
0467     // Separate spin up and down centroid positions
0468     float south_x_corr = smd_pos_corr[3];
0469     float south_y_corr = smd_pos_corr[2];
0470     float south_phi_corr = atan2(south_y_corr, south_x_corr);
0471     int yellowSpin = spinPatternYellow[bunchNum];
0472     if (yellowSpin == 1)
0473     {
0474       smd_hor_south_corrected_up->Fill(south_y_corr);
0475       smd_ver_south_corrected_up->Fill(south_x_corr);
0476       smd_south_phi_up_corrected->Fill(south_phi_corr);
0477       smd_south_phi_L_up_corrected->Fill(south_phi_corr);
0478       smd_south_phi_R_up_corrected->Fill(south_phi_corr);
0479     }
0480     if (yellowSpin == -1)
0481     {
0482       smd_hor_south_corrected_down->Fill(south_y_corr);
0483       smd_ver_south_corrected_down->Fill(south_x_corr);
0484       smd_south_phi_down_corrected->Fill(south_phi_corr);
0485       smd_south_phi_L_down_corrected->Fill(south_phi_corr);
0486       smd_south_phi_R_down_corrected->Fill(south_phi_corr);
0487     }
0488   }
0489
0490   // Now update the ADC values for the next event
0491   towerinfosZDC = findNode::getClass<TowerInfoContainer>(topNode, "TOWERS_ZDC");
0492   packetsZDC = findNode::getClass<CaloPacketContainerv1>(topNode, "ZDCPackets");
0493   GetAdcs();
0494   /* smdevent++; */
0495   /* CompSmdPos(); */
0496   /* std::cout << Form("Event %d end: smd_adc = {", evtctr); */
0497   /* for (int i=0; i<15; i++) { */
0498   /*     std::cout << smd_adc[i] << ", "; */
0499   /* } */
0500   /* std::cout << smd_adc[15] << "}" << std::endl; */
0501   /* std::cout << Form("Event %d end: smd_time = {", evtctr); */
0502   /* for (int i=0; i<15; i++) { */
0503   /*     std::cout << smd_time[i] << ", "; */
0504   /* } */
0505   /* std::cout << smd_time[15] << "}" << std::endl; */
0506   /* Gl1Packetv1* gl1 = findNode::getClass<Gl1Packetv1>(topNode, "GL1Packet"); */
0507   /* if (gl1) */
0508   /* { */
0509   /*   bunchNum = gl1->getBunchNumber(); */
0510   /*   bunchNum = (bunchNum + crossingShift)%NBUNCHES; */
0511   /*   /1* std::cout << "Got bunch number = " << bunchNum << std::endl; *1/ */
0512   /* } */
0513
0514   return Fun4AllReturnCodes::EVENT_OK;
0515 }
0516
0517 //____________________________________________________________________________..
0518 int SmdHistGen::ResetEvent(PHCompositeNode *topNode)
0519 {
0520     /* std::cout << "SmdHistGen::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl; */
0521     return Fun4AllReturnCodes::EVENT_OK;
0522 }
0523
0524 //____________________________________________________________________________..
0525 int SmdHistGen::EndRun(const int runnumber)
0526 {
0527   std::cout << "SmdHistGen::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0528   return Fun4AllReturnCodes::EVENT_OK;
0529 }
0530
0531 //____________________________________________________________________________..
0532 int SmdHistGen::End(PHCompositeNode *topNode)
0533 {
0534   std::cout << "SmdHistGen::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0535   std::cout << "Processed " << evtctr << " total events" << std::endl;
0536
0537   // Raw asymmetry
0538   smd_north_phi_up->Sumw2();
0539   smd_north_phi_down->Sumw2();
0540   smd_north_phi_sum->Add(smd_north_phi_up, smd_north_phi_down, 1.0, 1.0);
0541   smd_north_phi_diff->Add(smd_north_phi_up, smd_north_phi_down, 1.0, -1.0);
0542   smd_south_phi_up->Sumw2();
0543   smd_south_phi_down->Sumw2();
0544   smd_south_phi_sum->Add(smd_south_phi_up, smd_south_phi_down, 1.0, 1.0);
0545   smd_south_phi_diff->Add(smd_south_phi_up, smd_south_phi_down, 1.0, -1.0);
0546   smd_north_phi_up_corrected->Sumw2();
0547   smd_north_phi_down_corrected->Sumw2();
0548   smd_north_phi_sum_corrected->Add(smd_north_phi_up_corrected, smd_north_phi_down_corrected, 1.0, 1.0);
0549   smd_north_phi_diff_corrected->Add(smd_north_phi_up_corrected, smd_north_phi_down_corrected, 1.0, -1.0);
0550   smd_south_phi_up_corrected->Sumw2();
0551   smd_south_phi_down_corrected->Sumw2();
0552   smd_south_phi_sum_corrected->Add(smd_south_phi_up_corrected, smd_south_phi_down_corrected, 1.0, 1.0);
0553   smd_south_phi_diff_corrected->Add(smd_south_phi_up_corrected, smd_south_phi_down_corrected, 1.0, -1.0);
0554
0555   // Plot some histograms
0556   /*
0557   gStyle->SetOptStat(0);
0558   TCanvas* c1 = new TCanvas("c1", "c1", 800, 800);
0559   c1->cd();
0560   smd_hor_north->SetLineColor(kBlack);
0561   smd_hor_north_up->SetLineColor(kBlue);
0562   smd_hor_north_down->SetLineColor(kRed);
0563   smd_hor_north->Draw();
0564   smd_hor_north_up->Draw("same");
0565   smd_hor_north_down->Draw("same");
0566   TLegend* leg = new TLegend(0.7, 0.7, 0.9, 0.9);
0567   leg->AddEntry(smd_hor_north, "Total");
0568   leg->AddEntry(smd_hor_north_up, "Spin up");
0569   leg->AddEntry(smd_hor_north_down, "Spin down");
0570   leg->Draw();
0571   c1->Update();
0572   c1->SaveAs("plots/smd_hor_north_42200.png");
0573   delete leg;
0574
0575   smd_ver_north->SetLineColor(kBlack);
0576   smd_ver_north_up->SetLineColor(kBlue);
0577   smd_ver_north_down->SetLineColor(kRed);
0578   smd_ver_north->Draw();
0579   smd_ver_north_up->Draw("same");
0580   smd_ver_north_down->Draw("same");
0581   leg = new TLegend(0.7, 0.7, 0.9, 0.9);
0582   leg->AddEntry(smd_ver_north, "Total");
0583   leg->AddEntry(smd_ver_north_up, "Spin up");
0584   leg->AddEntry(smd_ver_north_down, "Spin down");
0585   leg->Draw();
0586   c1->Update();
0587   c1->SaveAs("plots/smd_ver_north_42200.png");
0588   delete leg;
0589
0590   smd_hor_north_neutron_hits->Draw();
0591   c1->Update();
0592   c1->SaveAs("plots/smd_hor_neutron_hits_42200.png");
0593   smd_ver_north_neutron_hits->Draw();
0594   c1->Update();
0595   c1->SaveAs("plots/smd_ver_neutron_hits_42200.png");
0596
0597   smd_north_waveforms->Draw("COLZ");
0598   gPad->SetLogz();
0599   c1->Update();
0600   c1->SaveAs("plots/smd_waveforms_42200.png");
0601   zdc_north_waveforms->Draw("COLZ");
0602   gPad->SetLogz();
0603   c1->Update();
0604   c1->SaveAs("plots/zdc_waveforms_42200.png");
0605   veto_north_waveforms->Draw("COLZ");
0606   gPad->SetLogz();
0607   c1->Update();
0608   c1->SaveAs("plots/veto_waveforms_42200.png");
0609   */
0610
0611   outfile->cd();
0612   outfile->Write();
0613   outfile->Close();
0614
0615   return Fun4AllReturnCodes::EVENT_OK;
0616 }
0617
0618 //____________________________________________________________________________..
0619 int SmdHistGen::Reset(PHCompositeNode *topNode)
0620 {
0621   std::cout << "SmdHistGen::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0622   return Fun4AllReturnCodes::EVENT_OK;
0623 }
0624
0625 //____________________________________________________________________________..
0626 void SmdHistGen::Print(const std::string &what) const
0627 {
0628    std::cout << "SmdHistGen::Print(const std::string &what) const Printing info for " << what << std::endl;
0629 }
0630
0631 int SmdHistGen::GetSpinPatterns()
0632 {
0633   // Get the spin patterns from the spin DB
0634
0635   //  0xffff is the qa_level from XingShiftCal //////
0636   unsigned int qa_level = 0xffff;
0637   SpinDBContent spin_cont;
0638   SpinDBOutput spin_out("phnxrc");
0639
0640   spin_out.StoreDBContent(runNum,runNum,qa_level);
0641   spin_out.GetDBContentStore(spin_cont,runNum);
0642
0643   // Get crossing shift
0644   crossingShift = spin_cont.GetCrossingShift();
0645   if (crossingShift == -999) crossingShift = 0;
0646   std::cout << "Crossing shift: " << crossingShift << std::endl;
0647
0648   std::cout << "Blue spin pattern: [";
0649   for (int i = 0; i < 120; i++)
0650   {
0651     spinPatternBlue[i] = spin_cont.GetSpinPatternBlue(i);
0652     std::cout << spinPatternBlue[i];
0653     if (i < 119)std::cout << ", ";
0654   }
0655   std::cout << "]" << std::endl;
0656
0657   std::cout << "Yellow spin pattern: [";
0658   for (int i = 0; i < 120; i++)
0659   {
0660     spinPatternYellow[i] = spin_cont.GetSpinPatternYellow(i);
0661     std::cout << spinPatternYellow[i];
0662     if (i < 119)std::cout << ", ";
0663   }
0664   std::cout << "]" << std::endl;
0665
0666   if (spinPatternYellow[0] == -999) {return 1;}
0667   return 0;
0668 }
0669
0670 void SmdHistGen::GetAdcs() // populate the adc arrays
0671 {
0672   if (mode == "dst") {GetAdcsDst();}
0673   else if (mode == "raw") {GetAdcsRaw();}
0674   else
0675   {
0676     std::cout << PHWHERE << ":: invalid mode=" << mode << " ... Exiting!" << std::endl;
0677     exit(1);
0678   }
0679 }
0680
0681 void SmdHistGen::GetAdcsDst()
0682 {
0683   if(!towerinfosZDC)
0684   {
0685     std::cout << PHWHERE << ":: No TOWERS_ZDC!" << std::endl; exit(1);
0686   }
0687
0688   int nchannels_zdc = towerinfosZDC->size();
0689   /* std::cout << "towerinfosZDC has " << nchannels_zdc << " channels" << std::endl; */
0690   for (int channel = 0; channel < nchannels_zdc;channel++)
0691   {
0692     // Channel mapping: ZDCS: 0-8; ZDCN: 9-15; SMDN: 16-31; SMDS: 32-47; Veto Counter N: 48 (front); Veto Counter N: 49 (back);  Veto Counter S: 50 (front); Veto Counter S: 51 (back)
0693     float zdc_e = towerinfosZDC->get_tower_at_channel(channel)->get_energy();
0694     float zdc_t = towerinfosZDC->get_tower_at_channel(channel)->get_time();
0695
0696     if (channel < 16) // ZDC
0697     {
0698       // ZDC Mapping:
0699       // even high gain, odd low gain
0700       // 0,1 S1; 2,3 S2; 4,5 S3; 6,7 Ssum
0701       // 8,9 N1; 10,11 N2; 12,13 N3; 14,15 Nsum
0702       zdc_adc[channel] = zdc_e;
0703       zdc_time[channel] = zdc_t;
0704       if (channel == 0) zdc1_south->Fill(zdc_e);
0705       if (channel == 2) zdc2_south->Fill(zdc_e);
0706       if (channel==0 || channel==2 || channel==4) zdc_south_waveforms->Fill(zdc_t, zdc_e);
0707       if (channel == 8) zdc1_north->Fill(zdc_e);
0708       if (channel == 10) zdc2_north->Fill(zdc_e);
0709       if (channel==8 || channel==10 || channel==12) zdc_north_waveforms->Fill(zdc_t, zdc_e);
0710     }
0711     if (channel >= 16 && channel < 48) // SMD
0712     {
0713       // SMD Mapping:
0714       // 0-7 North H1-8; 8-14 North V1-7; 15 North sum
0715       // 16-23 South H1-8; 24-30 South V1-7; 31 South sum
0716       int smd_channel = channel - 16;
0717       if (smd_channel < 15)
0718       {
0719     smd_north_signals[smd_channel]->Fill(zdc_e);
0720     smd_north_waveforms->Fill(zdc_t, zdc_e);
0721     smd_north_waveforms_zoomed->Fill(zdc_t, zdc_e);
0722     smd_north_channel_waveforms[smd_channel]->Fill(zdc_t, zdc_e);
0723       }
0724       if (smd_channel > 15)
0725       {
0726     smd_south_signals[smd_channel-16]->Fill(zdc_e);
0727     smd_south_waveforms->Fill(zdc_t, zdc_e);
0728     smd_south_waveforms_zoomed->Fill(zdc_t, zdc_e);
0729     smd_south_channel_waveforms[smd_channel-16]->Fill(zdc_t, zdc_e);
0730       }
0731       // Set the ADC to 0 if we're below threshold or out of time
0732       if (zdc_e < minSMDcut) zdc_e = 0;
0733       if (smd_channel < 16)
0734       {
0735     if ((zdc_t < smd_north_t_low) || (zdc_t > smd_north_t_high)) zdc_e = 0;
0736       }
0737       else
0738       {
0739     if ((zdc_t < smd_south_t_low) || (zdc_t > smd_south_t_high)) zdc_e = 0;
0740       }
0741       smd_adc[smd_channel] = zdc_e;
0742       smd_time[smd_channel] = zdc_t;
0743     }
0744     if (channel >= 48 && channel < 52) // Veto
0745     {
0746       // Veto Mapping:
0747       // 0 N front; 1 N back; 2 S front; 3 S back
0748       int veto_channel = channel - 48;
0749       veto_adc[veto_channel] = zdc_e;
0750       veto_time[veto_channel] = zdc_t;
0751       if (veto_channel == 0) vetofront_north->Fill(zdc_e);
0752       if (veto_channel == 1) vetoback_north->Fill(zdc_e);
0753       if (veto_channel < 2) veto_north_waveforms->Fill(zdc_t, zdc_e);
0754     }
0755   }
0756 }
0757
0758 void SmdHistGen::GetAdcsRaw()
0759 {
0760   if(!packetsZDC)
0761   {
0762     std::cout << PHWHERE << ":: No ZDCPackets!" << std::endl; exit(1);
0763   }
0764
0765   /* packetsZDC->identify(); */
0766   CaloPacket *packetZDC = packetsZDC->getPacketbyId(packet_smd);
0767   /* packetZDC->identify(); */
0768
0769   int Nchannels = std::min(128, packetZDC->iValue(0, "CHANNELS"));
0770   for (int channel = 0; channel < Nchannels; channel++)
0771   {
0772     // Channel mapping: ZDCS: 0-8; ZDCN: 9-15; SMDN: 48-63; SMDS: 112-127; Veto Counter N: 16 (front); Veto Counter N: 17 (back);  Veto Counter S: 80 (front); Veto Counter S: 81 (back)
0773     std::vector<float> resultFast = anaWaveformFast(packetZDC, channel);  // fast waveform fitting
0774     float zdc_e = resultFast.at(0);
0775     float zdc_t = resultFast.at(1);
0776     float zdc_ped = resultFast.at(2);
0777
0778     if (channel < 16) // ZDC
0779     {
0780       zdc_adc[channel] = zdc_e;
0781       zdc_time[channel] = zdc_t;
0782       if (channel == 8) zdc1_north->Fill(zdc_e);
0783       if (channel == 10) zdc2_north->Fill(zdc_e);
0784       if (channel==8 || channel==10 || channel==12) zdc_north_waveforms->Fill(zdc_t, zdc_e);
0785       continue;
0786     }
0787
0788     // North Veto counters
0789     if (channel == 16)
0790     {
0791       veto_adc[0] = zdc_e;
0792       veto_time[0] = zdc_t;
0793       vetofront_north->Fill(zdc_e);
0794       veto_north_waveforms->Fill(zdc_t, zdc_e);
0795       continue;
0796     }
0797     if (channel == 17)
0798     {
0799       veto_adc[1] = zdc_e;
0800       veto_time[1] = zdc_t;
0801       vetoback_north->Fill(zdc_e);
0802       veto_north_waveforms->Fill(zdc_t, zdc_e);
0803       continue;
0804     }
0805
0806     // North SMD
0807     if (channel >= 48 && channel < 63)
0808     {
0809       int smd_channel = channel - 48;
0810       if (smd_channel < 15)
0811       {
0812     smd_north_signal->Fill(zdc_e);
0813     smd_north_signals[smd_channel]->Fill(zdc_e);
0814     smd_north_waveforms->Fill(zdc_t, zdc_e);
0815     smd_north_waveforms_zoomed->Fill(zdc_t, zdc_e);
0816     smd_north_channel_waveforms[smd_channel]->Fill(zdc_t, zdc_e);
0817     smd_north_pedestal->Fill(zdc_ped);
0818       }
0819       // Set the ADC to 0 if we're below threshold or out of time
0820       if (zdc_e < minSMDcut) zdc_e = 0;
0821       if ((zdc_t < smd_north_t_low) || (zdc_t > smd_north_t_high)) zdc_e = 0;
0822       smd_adc[smd_channel] = zdc_e;
0823       smd_time[smd_channel] = zdc_t;
0824       continue;
0825     }
0826
0827     // South Veto counters
0828     if (channel == 80)
0829     {
0830       veto_adc[2] = zdc_e;
0831       veto_time[2] = zdc_t;
0832       vetofront_south->Fill(zdc_e);
0833       veto_south_waveforms->Fill(zdc_t, zdc_e);
0834       continue;
0835     }
0836     if (channel == 81)
0837     {
0838       veto_adc[3] = zdc_e;
0839       veto_time[3] = zdc_t;
0840       vetoback_south->Fill(zdc_e);
0841       veto_south_waveforms->Fill(zdc_t, zdc_e);
0842       continue;
0843     }
0844
0845     // South SMD
0846     if (channel >= 112 && channel < 127)
0847     {
0848       int smd_channel = channel - 112;
0849       if (smd_channel < 15)
0850       {
0851     smd_south_signal->Fill(zdc_e);
0852     smd_south_signals[smd_channel]->Fill(zdc_e);
0853     smd_south_waveforms->Fill(zdc_t, zdc_e);
0854     smd_south_waveforms_zoomed->Fill(zdc_t, zdc_e);
0855     smd_south_channel_waveforms[smd_channel]->Fill(zdc_t, zdc_e);
0856     smd_south_pedestal->Fill(zdc_ped);
0857       }
0858       // Set the ADC to 0 if we're below threshold or out of time
0859       if (zdc_e < minSMDcut) zdc_e = 0;
0860       if ((zdc_t < smd_south_t_low) || (zdc_t > smd_south_t_high)) zdc_e = 0;
0861       smd_adc[smd_channel+16] = zdc_e;
0862       smd_time[smd_channel+16] = zdc_t;
0863       continue;
0864     }
0865   }
0866 }
0867
0868 void SmdHistGen::CompSmdAdc() // mulitplying by relative gains
0869 {
0870   for (int i = 0; i < 15; i++) // last one is reserved for analogue sum
0871   {
0872     // multiply SMD channels with their gain factors
0873     // to get the absolute ADC values in the same units
0874     //rgains come from CompSmdAdc()
0875     smd_adc[i] = smd_adc[i] * smd_north_rgain[i]; // sout -> north for PHENIX -> sPHENIX
0876     smd_adc[i + 16] = smd_adc[i + 16] * smd_south_rgain[i]; // north -> south for PHENIX-> sPHENIX
0877
0878     /* for (int i=0; i<32; i++) { */
0879     /*   std::cout << "smd_adc[" << i << "] = " << smd_adc[i] << ", smd_time[" << i << "] = " << smd_time[i] << std::endl; */
0880     /* } */
0881     /* std::cout << std::endl; */
0882
0883   }
0884 }
0885
0886 void SmdHistGen::CompSmdPos() //computing position with weighted averages
0887 {
0888   float w_ave[4]; // 0 -> north hor; 1 -> noth vert; 2 -> south hor; 3 -> south vert.
0889   float weights[4] = {0};
0890   memset(weights, 0, sizeof(weights)); // memset float works only for 0
0891   float w_sum[4];
0892   memset(w_sum, 0, sizeof(w_sum));
0893
0894
0895   // these constants convert the SMD channel number into real dimensions (cm's)
0896   const float hor_scale = 2.0 * 11.0 / 10.5 * sin(PI/4); // from gsl_math.h
0897   const float ver_scale = 1.5 * 11.0 / 10.5;
0898   float hor_offset = (hor_scale * 8 / 2.0) * (7.0 / 8.0);
0899   float ver_offset = (ver_scale * 7 / 2.0) * (6.0 / 7.0);
0900
0901   for (int i = 0; i < 8; i++)
0902   {
0903     weights[0] += smd_adc[i]; //summing weights
0904     weights[2] += smd_adc[i + 16];
0905     w_sum[0] += (float)i * smd_adc[i]; //summing for the average
0906     w_sum[2] += ((float)i + 16.) * smd_adc[i + 16];
0907   }
0908   for (int i = 0; i < 7; i++)
0909   {
0910     weights[1] += smd_adc[i + 8];
0911     weights[3] += smd_adc[i + 24];
0912     w_sum[1] += ((float)i + 8.) * smd_adc[i + 8];
0913     w_sum[3] += ((float)i + 24.) * smd_adc[i + 24];
0914   }
0915
0916   if ( weights[0] > 0.0 )
0917   {
0918     w_ave[0] = w_sum[0] / weights[0]; //average = sum / sumn of weights...
0919     smd_pos[0] = hor_scale * w_ave[0] - hor_offset;
0920   }
0921   else
0922   {
0923     smd_pos[0] = 0;
0924   }
0925   if ( weights[1] > 0.0 )
0926   {
0927     w_ave[1] = w_sum[1] / weights[1];
0928     smd_pos[1] = ver_scale * (w_ave[1] - 8.0) - ver_offset;
0929   }
0930   else
0931   {
0932     smd_pos[1] = 0;
0933   }
0934
0935   if ( weights[2] > 0.0 )
0936   {
0937     w_ave[2] = w_sum[2] / weights[2];
0938     smd_pos[2] = hor_scale * (w_ave[2] - 16.0) - hor_offset;
0939   }
0940   else
0941   {
0942     smd_pos[2] = 0;
0943   }
0944
0945   if ( weights[3] > 0.0 )
0946   {
0947     w_ave[3] = w_sum[3] / weights[3];
0948     smd_pos[3] = ver_scale * (w_ave[3] - 24.0) - ver_offset;
0949   }
0950   else
0951   {
0952     smd_pos[3] = 0;
0953   }
0954 }
0955
0956 void SmdHistGen::CompSmdPosCorr()
0957 {
0958   /* float north_x_ctr = 0.351; */
0959   /* float north_y_ctr = 0.754; */
0960   /* float south_x_ctr = -0.334; */
0961   /* float south_y_ctr = -0.067; */
0962   float north_x_ctr = 0.407-0.031;
0963   float north_y_ctr = 0.790+0.149;
0964   float south_x_ctr = -0.393-0.176-0.062-0.019;
0965   float south_y_ctr = 0.082+0.037+0.014;
0966
0967   smd_pos_corr[1] = smd_pos[1] - north_x_ctr;
0968   smd_pos_corr[0] = smd_pos[0] - north_y_ctr;
0969   smd_pos_corr[3] = smd_pos[3] - south_x_ctr;
0970   smd_pos_corr[2] = smd_pos[2] - south_y_ctr;
0971 }
0972
0973 void SmdHistGen::CompSumSmd() //compute 'digital' sum
0974 {
0975   memset(smd_sum, 0, sizeof(smd_sum));
0976
0977   for (int i = 0; i < 8; i++)
0978   {
0979     smd_sum[0] += smd_adc[i]; // north horizontal
0980     smd_sum[2] += smd_adc[i + 16]; // south horizontal
0981   }
0982   for (int i = 0; i < 7; i++)
0983   {
0984     smd_sum[1] += smd_adc[i + 8]; // north vertical
0985     smd_sum[3] += smd_adc[i + 24]; // south vertical
0986   }
0987 }
0988
0989 void SmdHistGen::CountSMDHits()
0990 {
0991   // reset multiplicities
0992   n_hor_numhits = 0;
0993   n_ver_numhits = 0;
0994   s_hor_numhits = 0;
0995   s_ver_numhits = 0;
0996
0997   std::vector<int> n_hor(8, -999);
0998   std::vector<int> n_ver(7, -999);
0999   std::vector<int> s_hor(8, -999);
1000   std::vector<int> s_ver(7, -999);
1001   int n_hor_adjacent;
1002   int n_ver_adjacent;
1003   int s_hor_adjacent;
1004   int s_ver_adjacent;
1005
1006   // north
1007   for ( int i = 0; i < 8; i++)
1008   {
1009     n_hor[i] = 10*(i+1);
1010     if ( smd_adc[i] > minSMDcut )
1011     {
1012       // timing requirement
1013       if (smd_time[i]>=smd_north_t_low && smd_time[i]<=smd_north_t_high)
1014       {
1015         n_hor_numhits ++;
1016         smd_hor_north_total_hits->Fill(i);
1017     n_hor[i] = 0;
1018       }
1019     }
1020   }
1021   smd_hor_north_multiplicity->Fill(n_hor_numhits);
1022   for ( int i = 0; i < 7; i++)
1023   {
1024     n_ver[i] = 10*(i+1);
1025     if ( smd_adc[i + 8] > minSMDcut )
1026     {
1027       // timing requirement
1028       if (smd_time[i+8]>smd_north_t_low && smd_time[i+8]<=smd_north_t_high)
1029       {
1030         n_ver_numhits ++;
1031     smd_ver_north_total_hits->Fill(i);
1032     n_ver[i] = 0;
1033       }
1034     }
1035   }
1036   smd_ver_north_multiplicity->Fill(n_ver_numhits);
1037
1038   // south
1039   for ( int i = 0; i < 8; i++)
1040   {
1041     s_hor[i] = 10*(i+1);
1042     if ( smd_adc[i + 16] > minSMDcut )
1043     {
1044       if ( smd_time[i+16]>=smd_south_t_low && smd_time[i+16]<=smd_south_t_high )
1045       {
1046     s_hor_numhits++;
1047     smd_hor_south_total_hits->Fill(i);
1048     s_hor[i] = 0;
1049       }
1050     }
1051   }
1052   smd_hor_south_multiplicity->Fill(s_hor_numhits);
1053   for ( int i = 0; i < 7; i++)
1054   {
1055     s_ver[i] = 10*(i+1);
1056     if ( smd_adc[i + 24] > minSMDcut )
1057     {
1058       if ( smd_time[i+24]>=smd_south_t_low && smd_time[i+24]<=smd_south_t_high )
1059       {
1060     s_ver_numhits++;
1061     smd_ver_south_total_hits->Fill(i);
1062     s_ver[i] = 0;
1063       }
1064     }
1065   }
1066   smd_ver_south_multiplicity->Fill(s_ver_numhits);
1067
1068   n_hor_adjacent = CountAdjacentHits(n_hor);
1069   n_ver_adjacent = CountAdjacentHits(n_ver);
1070   s_hor_adjacent = CountAdjacentHits(s_hor);
1071   s_ver_adjacent = CountAdjacentHits(s_ver);
1072   // for testing
1073   if (false)
1074   {
1075     n_hor_numhits = n_hor_adjacent;
1076     n_ver_numhits = n_ver_adjacent;
1077     s_hor_numhits = s_hor_adjacent;
1078     s_ver_numhits = s_ver_adjacent;
1079   }
1080 }
1081
1082 bool SmdHistGen::NeutronSelection(std::string which, bool centerCorrected)
1083 {
1084   // Requirements:
1085   // veto ADC<200
1086   // ZDC2 ADC>20 (use high gain channels)
1087   // num SMD hits > 1
1088   // all hits in correct time window (north ZDC&Veto 5-8, north SMD 9-12)
1089   int frontveto, backveto, zdc1, zdc2, smdhitshor, smdhitsver;
1090   int frontveto_t, backveto_t, zdc1_t, zdc2_t, veto_t_low, veto_t_high;
1091   float smd_x, smd_y;
1092   if (which == "north") {
1093     frontveto = veto_adc[0];
1094     backveto = veto_adc[1];
1095     zdc1 = zdc_adc[8];
1096     zdc2 = zdc_adc[10];
1097     frontveto_t = veto_time[0];
1098     backveto_t = veto_time[1];
1099     zdc1_t = zdc_time[8];
1100     zdc2_t = zdc_time[10];
1101     smdhitshor = n_hor_numhits;
1102     smdhitsver = n_ver_numhits;
1103     veto_t_low = veto_north_t_low;
1104     veto_t_high = veto_north_t_high;
1105     smd_x = smd_pos[1];
1106     smd_y = smd_pos[0];
1107     if (centerCorrected)
1108     {
1109       smd_x = smd_pos_corr[1];
1110       smd_y = smd_pos_corr[0];
1111     }
1112   }
1113   else if (which == "south") {
1114     frontveto = veto_adc[2];
1115     backveto = veto_adc[3];
1116     zdc1 = zdc_adc[0];
1117     zdc2 = zdc_adc[2];
1118     frontveto_t = veto_time[2];
1119     backveto_t = veto_time[3];
1120     zdc1_t = zdc_time[0];
1121     zdc2_t = zdc_time[2];
1122     smdhitshor = s_hor_numhits;
1123     smdhitsver = s_ver_numhits;
1124     veto_t_low = veto_south_t_low;
1125     veto_t_high = veto_south_t_high;
1126     smd_x = smd_pos[3];
1127     smd_y = smd_pos[2];
1128     if (centerCorrected)
1129     {
1130       smd_x = smd_pos_corr[3];
1131       smd_y = smd_pos_corr[2];
1132     }
1133   }
1134   else {
1135     std::cout << "NeutronSelection: invalid string " << which << std::endl;
1136     return false;
1137   }
1138
1139   // timing requirement
1140   if (zdc1_t<zdc_t_low || zdc1_t>zdc_t_high) {return false;}
1141   if (zdc2_t<zdc_t_low || zdc2_t>zdc_t_high) {return false;}
1142   if (frontveto_t<veto_t_low || frontveto_t>veto_t_high) {return false;}
1143   if (backveto_t<veto_t_low || backveto_t>veto_t_high) {return false;}
1144   // ADC requirements
1145   if (frontveto > maxVetocut) {return false;}
1146   if (backveto > maxVetocut) {return false;}
1147   if (zdc1 < minZDC1cut) {return false;}
1148   if (zdc2 < minZDC2cut) {return false;}
1149   // SMD hit requirement
1150   if (smdhitshor < minSMDhits) {return false;}
1151   if (smdhitsver < minSMDhits) {return false;}
1152   // Radial position cut
1153   float r = sqrt(smd_x*smd_x + smd_y*smd_y);
1154   if (r < radius_low) {return false;}
1155   if (r > radius_high) {return false;}
1156   // passed all cuts
1157   return true;
1158 }
1159
1160 std::vector<float> SmdHistGen::anaWaveformFast(CaloPacket *p, const int channel)
1161 {
1162   std::vector<float> waveform;
1163   for (int s = 0; s < p->iValue(0, "SAMPLES"); s++)
1164   {
1165     waveform.push_back(p->iValue(s, channel));
1166     // fill waveform histograms
1167     if (channel >= 48 && channel < 63)
1168     {
1169       smd_north_waveforms_raw->Fill(s, p->iValue(s, channel));
1170     }
1171     if (channel >= 112 && channel < 127)
1172     {
1173       smd_south_waveforms_raw->Fill(s, p->iValue(s, channel));
1174     }
1175   }
1176   std::vector<std::vector<float>> multiple_wfs;
1177   multiple_wfs.push_back(waveform);
1178
1179   std::vector<std::vector<float>> fitresults_zdc;
1180   fitresults_zdc = WaveformProcessingFast->calo_processing_fast(multiple_wfs);
1181
1182   std::vector<float> result;
1183   result = fitresults_zdc.at(0);
1184   return result;
1185 }
1186
1187 int SmdHistGen::CountAdjacentHits(std::vector<int> channels)
1188 {
1189   // this function assumes each channel with a valid hit has value 0,
1190   // all other channels have *distinct* nonzero values
1191   std::vector<int> hit_indices;
1192   for (unsigned int i=0; i<channels.size(); i++)
1193   {
1194     if (channels[i] == 0) {hit_indices.push_back(i);}
1195   }
1196   int nhits = hit_indices.size();
1197   // if there is only 1 hit, let's say it is adjacent to itself
1198   if (nhits < 2) {return nhits;}
1199   int n_adjacent = 1;
1200   for (int i=0; i<(nhits-1); i++)
1201   {
1202     int delta = hit_indices[i+1] - hit_indices[i];
1203     if (delta == 1) {n_adjacent++;}
1204   }
1205   return n_adjacent;
1206 }