TemplateCreation/src/TemplateCreation.cc

0001
0002 #include "TemplateCreation.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <fun4all/SubsysReco.h>  // for SubsysReco
0006 #include <fun4all/PHTFileServer.h>
0007
0008
0009 #include <Event/Event.h>
0010 #include <Event/packet.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHIODataNode.h>    // for PHIODataNode
0014 #include <phool/PHNode.h>          // for PHNode
0015 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
0016 #include <phool/PHObject.h>        // for PHObject
0017 #include <phool/getClass.h>
0018
0019 #include <TProfile.h>
0020 #include <TMath.h>
0021 #include <TF1.h>
0022 #include <TH2.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025
0026
0027
0028
0029
0030
0031
0032 double TemplateCreation::rising_shape(double *val, double *par){
0033   double ped = par[0];
0034   double N_lan = par[1];
0035   double peak = par[2];
0036   double width = par[3];
0037   double x = val[0];
0038   double result = ped + N_lan*TMath::Landau(x, peak, width);
0039   return result;
0040 }
0041
0042 //____________________________________________________________________________..
0043 TemplateCreation::TemplateCreation(const std::string &name):
0044  SubsysReco(name)
0045  , s_outfilename("testout.root")
0046  , nsamples(31)
0047  , m_packet_low(21351)
0048  , m_packet_high(21352)
0049  , m_nchannels(192)
0050  , target_hm(5)
0051 {
0052 }
0053
0054 //____________________________________________________________________________..
0055 TemplateCreation::~TemplateCreation()
0056 {
0057 }
0058
0059 //____________________________________________________________________________..
0060 int TemplateCreation::Init(PHCompositeNode *topNode)
0061 {
0062   PHTFileServer::get().open(s_outfilename, "RECREATE");
0063   hp_all_particle_fine = new TProfile("hp_all_particle_fine","", nsamples*binscale, -0.5, nsamples - .5);
0064   h2_all_particle_fine = new TH2D("h2_all_particle_fine","", nsamples*binscale, -0.5, nsamples - .5, 100, -0.1, 0.5);
0065   return Fun4AllReturnCodes::EVENT_OK;
0066 }
0067
0068 //____________________________________________________________________________..
0069 int TemplateCreation::InitRun(PHCompositeNode *topNode)
0070 {
0071   return Fun4AllReturnCodes::EVENT_OK;
0072 }
0073
0074 //____________________________________________________________________________..
0075 int TemplateCreation::process_event(PHCompositeNode *topNode)
0076 {
0077   TH1F *h_waveform = new TH1F("h_waveform","", nsamples, -0.5, nsamples - 0.5);
0078   TF1 *f_rising_func = new TF1("f_rising_func",this,&TemplateCreation::rising_shape, 0, 31, 4,"TemplateCreation","rising_shape");
0079
0080   Event *_event = findNode::getClass<Event>(topNode, "PRDF");
0081   if (_event == nullptr)
0082     {
0083       std::cout << "CaloUnpackPRDF::Process_Event - Event not found" << std::endl;
0084       return -1;
0085     }
0086     if (_event->getEvtType() >= 8)  /// special event where we do not read out the calorimeters
0087     {
0088       return Fun4AllReturnCodes::DISCARDEVENT;
0089     }
0090     for (int pid = m_packet_low; pid <= m_packet_high; pid++)
0091     {
0092       Packet *packet = _event->getPacket(pid);
0093       if (packet)
0094     {
0095       int nchannels = packet->iValue(0, "CHANNELS");
0096       if (nchannels > m_nchannels) // packet is corrupted and reports too many channels
0097         {
0098           return Fun4AllReturnCodes::DISCARDEVENT;
0099         }
0100       for (int channel = 0; channel < nchannels; channel++)
0101         {
0102           float hm = 0;
0103           float sum = 0;
0104           int maxbin = 0;
0105           float shift = 0;
0106
0107           std::vector<float> waveform;
0108           waveform.reserve(nsamples);
0109           for (int samp = 0; samp < nsamples; samp++)
0110         {
0111           h_waveform->SetBinContent(samp + 1, packet->iValue(samp, channel));
0112         }
0113
0114           maxbin = h_waveform->GetMaximumBin();
0115           if (maxbin > 12){continue;} //Selection region really only needed for test beam data
0116           if (maxbin <= 5) {continue;}
0117           if (h_waveform->GetMaximum() < 2000){continue;}
0118
0119           float pedestal = 1500;
0120           if (maxbin > 5)
0121         {
0122           pedestal = 0.5 * (h_waveform->GetBinContent(maxbin - 6) + h_waveform->GetBinContent(maxbin - 5));
0123         }
0124
0125           for (int j = 1; j <= nsamples; j++)
0126         {
0127           if (h_waveform->GetBinContent(j) - pedestal  > 0)
0128             {
0129               sum += h_waveform->GetBinContent(j) - pedestal;
0130             }
0131         }
0132           for (int j = 0; j <= nsamples; j++)
0133         {
0134           if (((h_waveform->GetBinContent(j+1)) - pedestal) > 0)
0135             {
0136               h_waveform->SetBinContent(j + 1, ((h_waveform->GetBinContent(j+1)) - pedestal)/sum);
0137             }
0138           else
0139             {
0140               h_waveform->SetBinContent(j + 1,0);
0141             }
0142         }
0143
0144           if (h_waveform->GetBinContent(h_waveform->GetMaximumBin()) > 0.33){continue;}
0145
0146           bool pileup = false;
0147           for (int k = maxbin+1; k <= nsamples;k++)
0148         {
0149           if ( h_waveform->GetBinContent(k) > 1.05 *  h_waveform->GetBinContent(k-1))
0150             {
0151               pileup = true;
0152             }
0153         }
0154          if (pileup){continue;}
0155
0156           f_rising_func->SetParameters(0, 0.1, static_cast<double>(h_waveform->GetBinCenter(maxbin)), 1);//, p3, p4, p5);
0157           f_rising_func->FixParameter(0, 0);
0158           f_rising_func->SetParLimits(2,static_cast<double>(h_waveform->GetBinCenter(maxbin - 1)), static_cast<double>(h_waveform->GetBinCenter(maxbin + 1)));
0159           f_rising_func->SetParLimits(3, 0.5, 3);
0160           f_rising_func->SetRange(static_cast<float>(h_waveform->GetBinCenter(maxbin-3)), static_cast<float>(h_waveform->GetBinCenter(maxbin+1)));
0161           h_waveform->Fit("f_rising_func","NRQ");
0162           hm = f_rising_func->GetParameter(2) - 1.5*f_rising_func->GetParameter(3);
0163           shift = target_hm - hm;
0164
0165           bool skipwaveform = false;
0166           for (int i = h_waveform->GetMaximumBin() + 5; i < nsamples ;i++)
0167         {
0168           if (h_waveform->GetMaximum() < binscale*h_waveform->GetBinContent(i))
0169             {
0170               skipwaveform = true;
0171               continue;
0172             }
0173         }
0174           if (skipwaveform){continue;}
0175           for (int j = 0; j < nsamples; j++)
0176         {
0177           hp_all_particle_fine->Fill(j + shift, h_waveform->GetBinContent(j+1));
0178           h2_all_particle_fine->Fill(j + shift, h_waveform->GetBinContent(j+1));
0179         }
0180         }
0181     }
0182     }
0183
0184
0185     delete h_waveform;
0186
0187   return Fun4AllReturnCodes::EVENT_OK;
0188 }
0189
0190 //____________________________________________________________________________..
0191 int TemplateCreation::ResetEvent(PHCompositeNode *topNode)
0192 {
0193   return Fun4AllReturnCodes::EVENT_OK;
0194 }
0195
0196 //____________________________________________________________________________..
0197 int TemplateCreation::EndRun(const int runnumber)
0198 {
0199   return Fun4AllReturnCodes::EVENT_OK;
0200 }
0201
0202 //____________________________________________________________________________..
0203 int TemplateCreation::End(PHCompositeNode *topNode)
0204 {
0205   hp_all_particle_fine->Scale(1/hp_all_particle_fine->GetMaximum());
0206
0207
0208
0209   for (int i = 0 ; i < hp_all_particle_fine->GetNbinsX();i++)
0210     {
0211       if ( hp_all_particle_fine->GetBinContent(i+1) < 0.001)
0212     {
0213       hp_all_particle_fine->SetBinContent(i+1,0);
0214     }
0215     }
0216
0217
0218
0219   PHTFileServer::get().cd(s_outfilename);
0220   hp_all_particle_fine->SetName("h_template");
0221   hp_all_particle_fine->Write();
0222   h2_all_particle_fine->Write();
0223
0224
0225   return Fun4AllReturnCodes::EVENT_OK;
0226 }
0227
0228 //____________________________________________________________________________..
0229 int TemplateCreation::Reset(PHCompositeNode *topNode)
0230 {
0231   return Fun4AllReturnCodes::EVENT_OK;
0232 }
0233