Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hcal.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hcal.C
1 #include <iostream>
2 #include "hcal.h"
3 #include <algorithm>
4 #include <cmath>
5 #include <TMath.h>
6 #include <TGraph.h>
7 #include <TCanvas.h>
8 #include <TF1.h>
9 #include <TH1F.h>
10 #include <TH2F.h>
11 #include <TFile.h>
12 #include <TROOT.h>
13 #include <TStyle.h>
14 #include <Event/EventTypes.h>
15 #include <Event/packetConstants.h>
16 #include <Event/packet.h>
17 #include <Event/packet_hbd_fpgashort.h>
18 
19 using namespace std;
20 
21 //_________________________________________
22 hcal::hcal()
23 {
24  //Default values
25  display = false;
26  verbosity = false;
27  polarity_positive = true;
28  PEDESTAL = 2048;
29  plot_max = PEDESTAL + 100.;
30  plot_min = 0.;
31  fit_analysis = true;
32 }
33 
34 //________________________________________
35 void hcal::Initialize()
36 {
37  TString name;
38  cout << "Reading " << NCH << " Channels." << endl;
39  NSAMPLES = 24;
40  double max = 3000;
41  for(int ich=0; ich<NCH; ich++)
42  {
43  gpulse[ich] = new TGraph(NSAMPLES);
44  name = "gch"; name += ich;
45  gpulse[ich]->SetName(name);
46  gpulse[ich]->SetMarkerStyle(20);
47  gpulse[ich]->SetMarkerSize(0.4);
48 
49  name = "Signal_";
50  name+= ich;
51  h_PulsePeaks[ich] = new TH1F(name, name, (int)4*max, 0., max);
52  h_PulsePeaks[ich]->SetStats(kFALSE);
53 
54  //Pulse Int
55  name = "SignalInt_";
56  name+= ich;
57  double int_max = max*24./6. ;
58  h_PulseInt[ich] = new TH1F(name, name, (int)int_max, 0., int_max );
59 
60  //2D Integrated signal vs signal height
61  name = "Signal_Height_vs_Int_";
62  name += ich;
63  h2_Pulse_H_Int[ich] = new TH2F(name, name, (int) max, 0., max, (int)int_max, 0., int_max);
64 
65  //Fits
66  name = "SignalFit_";
67  name+= ich;
68  if(polarity_positive) fits[ich] = new TF1(name,hcal::SignalShape, 0., 24., 6);
69  else fits[ich] = new TF1(name,hcal::SignalShapeNegative, 0., 24., 6);
70 
71  //Fit analysis
72  if(fit_analysis)
73  {
74  name = "Fit_analyze_";
75  name += ich;
76  h2_fit_shape[ich] = new TH2F(name,name,48,0,24,2100,0,2100);
77  }
78 
79  //Timing resolution
80  name = "TResolution_";
81  name += ich;
82  h_tres[ich] = new TH1F(name, name, 24, 0, 24);
83 
84  //Pedestal
85  name = "Ped_Resolution_";
86  name += ich;
87  h_ped_res[ich] = new TH1F(name, name, 500,2000,2200);
88 
89  /*if(ich%2)
90  {
91  name = "HiLoGainRatio_";
92  name += (int)(ich/2);
93  h_gain[(int)ich/2] = new TH1F(name,name, 200, 0, 100);
94  name = "HiLo2D_";
95  name += (int)(ich/2);
96  h_hilo[(int)ich/2] = new TH2F(name,name, 500, 0, 3000, 500, 0, 300);
97  }*/
98  }
99 
100  name = "Info";
101  h_info = new TH1F(name, name, 10, 0, 10);
102  //1st bin = Runnumber
103  //2nd bin = Number of events ran
104  //3rd bin = Start time
105  //4th bin = End time
106  // More to add later
107 
108  name = "Saturation_rates_vs_ch";
109  h_saturation_rates = new TH1F(name, name, NCH, 0, NCH);
110 
111  cout << "Reading below HBD channels: " << endl;
112  for(int ich=0; ich<NCH; ich++) cout << channel_index[ich] << ", ";
113  cout << endl;
114 
115  //Save the trigger rate
116  rate = new TGraph();
117  rate->SetNameTitle("Rate","Trigger Rate");
118 }
119 
120 //____________________________________
121 int hcal::evLoop(int nevts=0)
122 {
123  Initialize();
124  runnumber = 0;
125  nevents = 0;
126  int OK;
127  if(verbosity) cout << "hcal::evLoop" << endl;
128  //Event iterator
129  cout << "Reading file " << prdfName << endl;
130  it = new fileEventiterator(prdfName, OK);
131  if (OK)
132  {
133  cout << "Couldn't open input file " << prdfName << endl;
134  delete(it);
135  exit(1);
136  }
137 
138  double time_begin = 0.;
139  double time_end = 0.;
140  while(GetNextEvent())
141  {
142  if(verbosity) cout << "Event " << nevents << endl;
143  if(display) Display();
144  if(nevts>0&&nevents==nevts) break;
145  if(nevents%100==0)
146  {
147  time_end = (double) evt->getTime();
148  if(nevents>100) rate->SetPoint((nevents/100)-2, nevents, 100./(time_end-time_begin));
149  time_begin = time_end;
150  }
151  }
152  cout << "Total number of events processed " << nevents << endl;
153  h_info->SetBinContent(2, nevents);
154 
155  cout << "Saturation summary " << endl;
156  for( int ich=0; ich<NCH; ich++)
157  {
158  cout << "Channel " << ich << " => " << 100*saturation_map[ich]/nevents << " % " <<endl;
159  h_saturation_rates->Fill( h_saturation_rates->FindBin(ich), 100*saturation_map[ich]/nevents );
160  }
161 
162 
163  return 0;
164 }
165 
166 //_____________________________________
167 int hcal::Clear()
168 {
169  if(verbosity) cout << "hcal::Clear()" << endl;
170  for (int ich=0; ich<NCH; ich++)
171  {
172  gpulse[ich]->Set(0);
173  //if(fits[ich]) delete fits[ich];
174  }
175  return 0;
176 }
177 
178 //______________________________________
179 int hcal::GetNextEvent()
180 {
181  if(verbosity) cout << "hcal::GetNextEvent " << nevents << endl;
182  nevents++;
183  //Reset the graphs
184  Clear();
185  //Get the event
186  evt=it->getNextEvent();
187  if(!evt) return 0;
188  if(evt->getEvtType()==9)
189  {
190  cout << "Event " << nevents << endl;
191  cout << "Start Time " << evt->getTime() << endl;
192  cout << "Start Date " << evt->getDate() << endl;
193  //Get Runnumber
194  h_info->SetBinContent(1, evt->getRunNumber() );
195  //Get Start Time
196  h_info->SetBinContent(3, evt->getTime() );
197  }
198  else if(evt->getEvtType()==12)
199  {
200  //Get End Time
201  h_info->SetBinContent(4, evt->getTime() );
202  } else {
203  // Get sPHENIX Packet
204  Packet_hbd_fpgashort *spacket = dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101) );
205  if ( spacket )
206  {
207  spacket->setNumSamples( NSAMPLES );
208  if(nevents%500==0) cout <<"RUN "<<runnumber<<" EVENT "<< nevents <<endl;
209  for(int ich=0; ich<NCH; ich++)
210  {
211  int hbd_channel = channel_index[ich];
212  for (int isamp=0; isamp<NSAMPLES; isamp++)
213  {
214  //Short_t val = spacket->iValue(feech[ich],isamp);
215  //gpulse[ich]->SetPoint(isamp,(Double_t)isamp,(Double_t)val);
216  Short_t val = spacket->iValue(hbd_channel,isamp);
217  gpulse[ich]->SetPoint(isamp,(Double_t)isamp,(Double_t)val);
218  }
219 
220  if(!polarity_positive && TMath::MinElement(24,gpulse[ich]->GetY())==0)
221  {
222  //display = true;
223  saturation_map[ich]++;
224  if(verbosity) cout << "Saturated channel " << ich << endl;
225  }
226  else
227  {
228  //display = false;
229  fitShape(ich);
230  collect(ich);
231  }
232  }
233  delete spacket;
234  } else {
235  if( nevents>1) cout << "Event: " << nevents << " Packet not found " << endl;
236  }
237  }
238  return 1;
239 }
240 
241 //____________________________________________
242 void hcal::fitShape(int ich)
243 {
244  int peakPos = 0.;
245  double peakval;
246  if(polarity_positive) peakval = 0.;
247  else peakval = 99999.;
248  double pedestal = gpulse[ich]->GetY()[0]; //(double) PEDESTAL;
249  double risetime = 4;
250  for(int iSample=0; iSample<NSAMPLES; iSample++)
251  {
252 
253  if(!polarity_positive && gpulse[ich]->GetY()[iSample]<peakval){
254  //if(gpulse[ich]->GetY()[iSample]>peakval){
255  peakval = gpulse[ich]->GetY()[iSample];
256  peakPos = iSample;
257  } else if(polarity_positive && gpulse[ich]->GetY()[iSample]>peakval)
258  {
259  peakval = gpulse[ich]->GetY()[iSample];
260  peakPos = iSample;
261  }
262  }
263 
264  //peakval = pedestal - peakval;
265  if(polarity_positive) peakval -= pedestal;
266  else peakval = pedestal - peakval;
267  if(peakval<0.) peakval = 0.;
268  double par[6];
269  par[0] = peakval; // /3.;
270  par[1] = peakPos - risetime;
271  if(par[1]<0.) par[1] = 0.;
272  par[2] = 4.;
273  par[3] = 1.5;
274  par[4] = pedestal;
275  par[5] = 0;
276  fits[ich]->SetParameters(par);
277  fits[ich]->SetParLimits(0,peakval*0.98,peakval*1.02);
278  fits[ich]->SetParLimits(1,0,24);
279  fits[ich]->SetParLimits(2, 2, 4.);
280  fits[ich]->SetParLimits(3,1.,2.);
281  fits[ich]->SetParLimits(4, pedestal-30, pedestal+30);
282  fits[ich]->SetParLimits(5, 1, -1);
283  gpulse[ich]->Fit( fits[ich], "MQR", "goff", 0., (double) NSAMPLES);
284  //gpulse[ich]->Fit( fits[ich], "M0RQ", "goff", 0., (double) NSAMPLES);
285  //cout << "CH " << ich << " sig: " << peakval << " Par 2: " << fits[ich]->GetParameter(2) << endl;
286  //fits[ich]->Print();
287 }
288 
289 //__________________________________________
290 int hcal::collect(int ich)
291 {
292  if(nevents<10) return 0; //Leave first 10 events
293  double min = fits[ich]->GetMinimum();
294  double minx = fits[ich]->GetMinimumX();
295  //double ped_val = fits[ich]->GetParameter(4)+minx*fits[ich]->GetParameter(5);
296  //double peakvalue = ped_val - min;
297 
298  double max = fits[ich]->GetMaximum();
299  double maxx = fits[ich]->GetMaximumX();
300  double ped_val =0;
301  if(polarity_positive) ped_val = fits[ich]->GetParameter(4)+maxx*fits[ich]->GetParameter(5);
302  else ped_val = fits[ich]->GetParameter(4)+minx*fits[ich]->GetParameter(5);
303 
304  double peakvalue = 0;
305  if(polarity_positive) peakvalue = max - ped_val;
306  else peakvalue = ped_val - min;
307 
308  double integral = 24*ped_val - fits[ich]->Integral(0,24);
309  if(verbosity) cout << "Channel " << ich << " " << h_PulsePeaks[ich]->GetName() << " PED: " << ped_val << " Signal " << peakvalue << " Integrated signal: " << integral << endl;
310 
311  //cout << "peakpos " << fits[ich]->GetMinimumX() << endl;
312  //cout << "Maximum " << fits[ich]->GetMaximum() << endl;
313  h_PulsePeaks[ich]->Fill( peakvalue );
314  h_PulseInt[ich]->Fill( integral );
315  h2_Pulse_H_Int[ich]->Fill( peakvalue, integral );
316  h_tres[ich]->Fill( (polarity_positive)?maxx:minx );
317  h_ped_res[ich]->Fill( ped_val );
318  /*if(ich%2)
319  {
320  double hi_peak = fits[ich-1]->GetParameter(4)+fits[ich-1]->GetMinimumX()*fits[ich-1]->GetParameter(5) - fits[ich-1]->GetMinimum();
321  double lo_peak = fits[ich-1]->GetMaximum() - fits[ich-1]->GetParameter(4)+fits[ich-1]->GetMaximumX()*fits[ich-1]->GetParameter(5);
322  if(polarity_positive)
323  {
324  h_gain[(int)ich/2]->Fill( peakvalue/lo_peak );
325  h_hilo[(int)ich/2]->Fill( peakvalue, lo_peak );
326  } else {
327  h_gain[(int)ich/2]->Fill( hi_peak/peakvalue );
328  h_hilo[(int)ich/2]->Fill( hi_peak, peakvalue );
329  }
330  //cout << "Channel " << ich << " " << (int)ich/2 << " " << hi_peak/peakvalue << endl;
331  }*/
332 
333  return 0;
334 }
335 
336 //_________________________________________
337 void hcal::Display()
338 {
339  if(nevents==1) return;
340  if(gpulse[0]->GetN()==0) return;
341 
342  if(fit_analysis)
343  {
344  for(int ich=0; ich<NCH; ich++)
345  {
346  if(verbosity) fits[ich]->Print();
347  fits[ich]->SetParameter(1,4.0);
348  for(float x=0; x<48; x=x+0.5)
349  {
350  h2_fit_shape[ich]->Fill( x, fits[ich]->Eval(x) );
351  }
352  }
353  }
354 
355  if(!c)
356  {
357  int nx_c = 4;
358  int ny_c = ceil( NCH/4 );
359  c = new TCanvas("sphenix","sphenix display", 300*nx_c, 200*ny_c );
360  c->Divide(nx_c,ny_c);
361  //int ncd[]={13,9,5,1,14,10,6,2,15,11,7,3,16,12,8,4};
362  for(int ich=0; ich<NCH; ich++)
363  {
364  //c->cd(ncd[ich]); //ich+1);
365  c->cd(ich+1);
366  gpulse[ich]->SetMaximum( plot_max ); //PEDESTAL+50 );
367  gpulse[ich]->SetMinimum( plot_min ); //PEDESTAL-2000 );
368  gpulse[ich]->Draw("ACP");
369  }
370 
371  //cout << "Double click on the canvas to go to next event." << endl;
372  cout << "Kill the process or 'killall root.exe' to stop the canvas" << endl;
373  }
374 
375  for(int ich=0; ich<NCH; ich++) c->cd(ich+1)->Modified();
376  //c->WaitPrimitive();
377  c->Modified();
378  c->Update();
379  //c->Print("print.gif+");
380 }
381 
382 //_____________________________________
383 double hcal::SignalShape(double *x, double *par)
384 {
385  double pedestal = par[4] + x[0]*par[5];
386  if( x[0]<par[1]) return pedestal;
387  //double signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
388  double signal = par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
389  return pedestal + signal ;
390 }
391 
392 //______________________________________
393 double hcal::SignalShapeNegative(double *x, double *par)
394 {
395  double pedestal = par[4] + x[0]*par[5];
396  if( x[0]<par[1]) return pedestal;
397  double signal = par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
398  return pedestal-signal ;
399 }
400 
401 
402 //________________________________________
403 void hcal::Save(char *filename)
404 {
405  fout = TFile::Open( filename, "RECREATE" );
406  for(int ich=0; ich<NCH; ich++)
407  {
408  h_PulsePeaks[ich]->Write();
409  h_tres[ich]->Write();
410  h_ped_res[ich]->Write();
411  h2_fit_shape[ich]->Write();
412  //h_hilo[ich/2]->Write();
413  //if(ich%2) h_gain[(int)ich/2]->Write();
414  h_PulseInt[ich]->Write();
415  h2_Pulse_H_Int[ich]->Write();
416  }
417  h_saturation_rates->Write();
418  rate->Write();
419  h_info->Write();
420  fout->Close();
421 }
422