HCAL/HighBayCode/hcal.C

0001 #include <iostream>
0002 #include "hcal.h"
0003 #include <algorithm>
0004 #include <cmath>
0005 #include <TMath.h>
0006 #include <TGraph.h>
0007 #include <TCanvas.h>
0008 #include <TF1.h>
0009 #include <TH1F.h>
0010 #include <TH2F.h>
0011 #include <TFile.h>
0012 #include <TROOT.h>
0013 #include <TStyle.h>
0014 #include <Event/EventTypes.h>
0015 #include <Event/packetConstants.h>
0016 #include <Event/packet.h>
0017 #include <Event/packet_hbd_fpgashort.h>
0018
0019 using namespace std;
0020
0021 //_________________________________________
0022 hcal::hcal()
0023 {
0024  //Default values
0025  display = false;
0026  verbosity = false;
0027  polarity_positive = true;
0028  PEDESTAL = 2048;
0029  plot_max = PEDESTAL + 100.;
0030  plot_min = 0.;
0031  fit_analysis = true;
0032 }
0033
0034 //________________________________________
0035 void hcal::Initialize()
0036 {
0037   TString name;
0038   cout << "Reading " << NCH << " Channels." << endl;
0039   NSAMPLES = 24;
0040   double max = 3000;
0041   for(int ich=0; ich<NCH; ich++)
0042     {
0043       gpulse[ich] = new TGraph(NSAMPLES);
0044       name = "gch"; name += ich;
0045       gpulse[ich]->SetName(name);
0046       gpulse[ich]->SetMarkerStyle(20);
0047       gpulse[ich]->SetMarkerSize(0.4);
0048
0049       name = "Signal_";
0050       name+= ich;
0051       h_PulsePeaks[ich] = new TH1F(name, name, (int)4*max, 0., max);
0052       h_PulsePeaks[ich]->SetStats(kFALSE);
0053
0054       //Pulse Int
0055       name = "SignalInt_";
0056       name+= ich;
0057       double int_max = max*24./6. ;
0058       h_PulseInt[ich] = new TH1F(name, name, (int)int_max, 0., int_max );
0059
0060       //2D Integrated signal vs signal height
0061       name = "Signal_Height_vs_Int_";
0062       name += ich;
0063       h2_Pulse_H_Int[ich] = new TH2F(name, name, (int) max, 0., max, (int)int_max, 0., int_max);
0064
0065       //Fits
0066       name = "SignalFit_";
0067       name+= ich;
0068       if(polarity_positive) fits[ich] = new TF1(name,hcal::SignalShape, 0., 24., 6);
0069       else fits[ich] = new TF1(name,hcal::SignalShapeNegative, 0., 24., 6);
0070
0071       //Fit analysis
0072       if(fit_analysis)
0073       {
0074        name = "Fit_analyze_";
0075        name += ich;
0076        h2_fit_shape[ich] = new TH2F(name,name,48,0,24,2100,0,2100);
0077       }
0078
0079       //Timing resolution
0080       name = "TResolution_";
0081       name += ich;
0082       h_tres[ich] = new TH1F(name, name, 24, 0, 24);
0083
0084       //Pedestal
0085       name = "Ped_Resolution_";
0086       name += ich;
0087       h_ped_res[ich] = new TH1F(name, name, 500,2000,2200);
0088
0089       /*if(ich%2)
0090       {
0091        name = "HiLoGainRatio_";
0092        name += (int)(ich/2);
0093        h_gain[(int)ich/2] = new TH1F(name,name, 200, 0, 100);
0094        name = "HiLo2D_";
0095        name += (int)(ich/2);
0096        h_hilo[(int)ich/2] = new TH2F(name,name, 500, 0, 3000, 500, 0, 300);
0097       }*/
0098     }
0099
0100   name = "Info";
0101   h_info = new TH1F(name, name, 10, 0, 10);
0102   //1st bin = Runnumber
0103   //2nd bin = Number of events ran
0104   //3rd bin = Start time
0105   //4th bin = End time
0106   // More to add later
0107
0108   name = "Saturation_rates_vs_ch";
0109   h_saturation_rates = new TH1F(name, name, NCH, 0, NCH);
0110
0111   cout << "Reading below HBD channels:  " << endl;
0112   for(int ich=0; ich<NCH; ich++) cout << channel_index[ich] << ", ";
0113   cout << endl;
0114
0115   //Save the trigger rate
0116   rate = new TGraph();
0117   rate->SetNameTitle("Rate","Trigger Rate");
0118 }
0119
0120 //____________________________________
0121 int hcal::evLoop(int nevts=0)
0122 {
0123  Initialize();
0124  runnumber = 0;
0125  nevents = 0;
0126  int OK;
0127  if(verbosity) cout << "hcal::evLoop" << endl;
0128  //Event iterator
0129  cout << "Reading file " << prdfName << endl;
0130  it =  new fileEventiterator(prdfName, OK);
0131  if (OK)
0132  {
0133   cout << "Couldn't open input file " << prdfName << endl;
0134   delete(it);
0135   exit(1);
0136  }
0137
0138  double time_begin = 0.;
0139  double time_end = 0.;
0140  while(GetNextEvent())
0141  {
0142   if(verbosity) cout << "Event " << nevents << endl;
0143   if(display) Display();
0144   if(nevts>0&&nevents==nevts) break;
0145   if(nevents%100==0)
0146   {
0147    time_end = (double) evt->getTime();
0148    if(nevents>100) rate->SetPoint((nevents/100)-2, nevents, 100./(time_end-time_begin));
0149    time_begin = time_end;
0150   }
0151  }
0152  cout << "Total number of events processed " << nevents << endl;
0153  h_info->SetBinContent(2, nevents);
0154
0155  cout << "Saturation summary " << endl;
0156  for( int ich=0; ich<NCH; ich++)
0157  {
0158   cout << "Channel " << ich << " => " << 100*saturation_map[ich]/nevents << " % " <<endl;
0159   h_saturation_rates->Fill( h_saturation_rates->FindBin(ich), 100*saturation_map[ich]/nevents );
0160  }
0161
0162
0163  return 0;
0164 }
0165
0166 //_____________________________________
0167 int hcal::Clear()
0168 {
0169  if(verbosity) cout << "hcal::Clear()" << endl;
0170   for (int ich=0; ich<NCH; ich++)
0171   {
0172     gpulse[ich]->Set(0);
0173     //if(fits[ich]) delete fits[ich];
0174   }
0175  return 0;
0176 }
0177
0178 //______________________________________
0179 int hcal::GetNextEvent()
0180 {
0181   if(verbosity) cout << "hcal::GetNextEvent " << nevents << endl;
0182   nevents++;
0183   //Reset the graphs
0184   Clear();
0185   //Get the event
0186   evt=it->getNextEvent();
0187   if(!evt) return 0;
0188   if(evt->getEvtType()==9)
0189   {
0190     cout << "Event " << nevents << endl;
0191     cout << "Start Time " << evt->getTime() << endl;
0192     cout << "Start Date " << evt->getDate() << endl;
0193     //Get Runnumber
0194     h_info->SetBinContent(1, evt->getRunNumber() );
0195     //Get Start Time
0196     h_info->SetBinContent(3, evt->getTime() );
0197   }
0198   else if(evt->getEvtType()==12)
0199   {
0200     //Get End Time
0201     h_info->SetBinContent(4, evt->getTime() );
0202   } else {
0203    // Get sPHENIX Packet
0204    Packet_hbd_fpgashort *spacket = dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101) );
0205    if ( spacket )
0206    {
0207      spacket->setNumSamples( NSAMPLES );
0208      if(nevents%500==0)  cout <<"RUN  "<<runnumber<<"  EVENT "<< nevents <<endl;
0209      for(int ich=0; ich<NCH; ich++)
0210      {
0211       int hbd_channel = channel_index[ich];
0212       for (int isamp=0; isamp<NSAMPLES; isamp++)
0213       {
0214        //Short_t val = spacket->iValue(feech[ich],isamp);
0215        //gpulse[ich]->SetPoint(isamp,(Double_t)isamp,(Double_t)val);
0216         Short_t val = spacket->iValue(hbd_channel,isamp);
0217         gpulse[ich]->SetPoint(isamp,(Double_t)isamp,(Double_t)val);
0218       }
0219
0220        if(!polarity_positive && TMath::MinElement(24,gpulse[ich]->GetY())==0)
0221        {
0222          //display = true;
0223          saturation_map[ich]++;
0224          if(verbosity) cout << "Saturated channel " << ich << endl;
0225        }
0226        else
0227        {
0228         //display = false;
0229         fitShape(ich);
0230         collect(ich);
0231        }
0232       }
0233        delete spacket;
0234     } else {
0235      if( nevents>1) cout << "Event: " << nevents << " Packet not found " << endl;
0236    }
0237   }
0238  return 1;
0239 }
0240
0241 //____________________________________________
0242 void hcal::fitShape(int ich)
0243 {
0244  int peakPos = 0.;
0245  double peakval;
0246  if(polarity_positive) peakval = 0.;
0247  else peakval = 99999.;
0248  double pedestal = gpulse[ich]->GetY()[0]; //(double) PEDESTAL;
0249  double risetime = 4;
0250  for(int iSample=0; iSample<NSAMPLES; iSample++)
0251  {
0252
0253    if(!polarity_positive && gpulse[ich]->GetY()[iSample]<peakval){
0254    //if(gpulse[ich]->GetY()[iSample]>peakval){
0255      peakval = gpulse[ich]->GetY()[iSample];
0256      peakPos = iSample;
0257    } else if(polarity_positive && gpulse[ich]->GetY()[iSample]>peakval)
0258    {
0259      peakval = gpulse[ich]->GetY()[iSample];
0260      peakPos = iSample;
0261    }
0262  }
0263
0264  //peakval = pedestal - peakval;
0265  if(polarity_positive) peakval -= pedestal;
0266  else peakval = pedestal - peakval;
0267  if(peakval<0.) peakval = 0.;
0268  double par[6];
0269  par[0] = peakval; // /3.;
0270  par[1] = peakPos - risetime;
0271  if(par[1]<0.) par[1] = 0.;
0272  par[2] = 4.;
0273  par[3] = 1.5;
0274  par[4] = pedestal;
0275  par[5] = 0;
0276  fits[ich]->SetParameters(par);
0277  fits[ich]->SetParLimits(0,peakval*0.98,peakval*1.02);
0278  fits[ich]->SetParLimits(1,0,24);
0279  fits[ich]->SetParLimits(2, 2, 4.);
0280  fits[ich]->SetParLimits(3,1.,2.);
0281  fits[ich]->SetParLimits(4, pedestal-30, pedestal+30);
0282  fits[ich]->SetParLimits(5, 1, -1);
0283  gpulse[ich]->Fit( fits[ich], "MQR", "goff", 0., (double) NSAMPLES);
0284  //gpulse[ich]->Fit( fits[ich], "M0RQ", "goff", 0., (double) NSAMPLES);
0285  //cout << "CH " << ich << " sig: " << peakval << " Par 2: " << fits[ich]->GetParameter(2) << endl;
0286  //fits[ich]->Print();
0287 }
0288
0289 //__________________________________________
0290 int hcal::collect(int ich)
0291 {
0292  if(nevents<10) return 0; //Leave first 10 events
0293  double min = fits[ich]->GetMinimum();
0294  double minx = fits[ich]->GetMinimumX();
0295  //double ped_val = fits[ich]->GetParameter(4)+minx*fits[ich]->GetParameter(5);
0296  //double peakvalue = ped_val - min;
0297
0298  double max = fits[ich]->GetMaximum();
0299  double maxx = fits[ich]->GetMaximumX();
0300  double ped_val =0;
0301  if(polarity_positive) ped_val = fits[ich]->GetParameter(4)+maxx*fits[ich]->GetParameter(5);
0302  else ped_val = fits[ich]->GetParameter(4)+minx*fits[ich]->GetParameter(5);
0303
0304  double peakvalue = 0;
0305  if(polarity_positive) peakvalue = max  - ped_val;
0306  else peakvalue = ped_val - min;
0307
0308  double integral = 24*ped_val - fits[ich]->Integral(0,24);
0309  if(verbosity) cout << "Channel " << ich << " " << h_PulsePeaks[ich]->GetName() << " PED: " << ped_val << " Signal " << peakvalue << " Integrated signal: " << integral << endl;
0310
0311  //cout << "peakpos " << fits[ich]->GetMinimumX() << endl;
0312  //cout << "Maximum " << fits[ich]->GetMaximum() << endl;
0313  h_PulsePeaks[ich]->Fill( peakvalue );
0314  h_PulseInt[ich]->Fill( integral );
0315  h2_Pulse_H_Int[ich]->Fill( peakvalue, integral );
0316  h_tres[ich]->Fill( (polarity_positive)?maxx:minx );
0317  h_ped_res[ich]->Fill( ped_val );
0318  /*if(ich%2)
0319  {
0320   double hi_peak = fits[ich-1]->GetParameter(4)+fits[ich-1]->GetMinimumX()*fits[ich-1]->GetParameter(5) - fits[ich-1]->GetMinimum();
0321   double lo_peak = fits[ich-1]->GetMaximum() - fits[ich-1]->GetParameter(4)+fits[ich-1]->GetMaximumX()*fits[ich-1]->GetParameter(5);
0322   if(polarity_positive)
0323   {
0324    h_gain[(int)ich/2]->Fill( peakvalue/lo_peak );
0325    h_hilo[(int)ich/2]->Fill( peakvalue, lo_peak );
0326   } else {
0327    h_gain[(int)ich/2]->Fill( hi_peak/peakvalue );
0328    h_hilo[(int)ich/2]->Fill( hi_peak, peakvalue );
0329   }
0330   //cout << "Channel " << ich << " " << (int)ich/2 << "  " << hi_peak/peakvalue << endl;
0331  }*/
0332
0333  return 0;
0334 }
0335
0336 //_________________________________________
0337 void hcal::Display()
0338 {
0339  if(nevents==1) return;
0340  if(gpulse[0]->GetN()==0) return;
0341
0342  if(fit_analysis)
0343  {
0344   for(int ich=0; ich<NCH; ich++)
0345   {
0346    if(verbosity) fits[ich]->Print();
0347    fits[ich]->SetParameter(1,4.0);
0348    for(float x=0; x<48; x=x+0.5)
0349    {
0350     h2_fit_shape[ich]->Fill( x, fits[ich]->Eval(x) );
0351    }
0352   }
0353  }
0354
0355  if(!c)
0356  {
0357   int nx_c = 4;
0358   int ny_c = ceil( NCH/4 );
0359   c = new TCanvas("sphenix","sphenix display", 300*nx_c, 200*ny_c );
0360   c->Divide(nx_c,ny_c);
0361   //int ncd[]={13,9,5,1,14,10,6,2,15,11,7,3,16,12,8,4};
0362   for(int ich=0; ich<NCH; ich++)
0363   {
0364    //c->cd(ncd[ich]); //ich+1);
0365    c->cd(ich+1);
0366    gpulse[ich]->SetMaximum( plot_max ); //PEDESTAL+50 );
0367    gpulse[ich]->SetMinimum( plot_min ); //PEDESTAL-2000 );
0368    gpulse[ich]->Draw("ACP");
0369   }
0370
0371   //cout << "Double click on the canvas to go to next event." << endl;
0372   cout << "Kill the process or 'killall root.exe' to stop the canvas" << endl;
0373  }
0374
0375  for(int ich=0; ich<NCH; ich++) c->cd(ich+1)->Modified();
0376  //c->WaitPrimitive();
0377  c->Modified();
0378  c->Update();
0379  //c->Print("print.gif+");
0380 }
0381
0382 //_____________________________________
0383 double hcal::SignalShape(double *x, double *par)
0384 {
0385  double pedestal = par[4] + x[0]*par[5];
0386  if( x[0]<par[1]) return pedestal;
0387  //double  signal = (-1)*par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0388  double  signal = par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0389  return pedestal + signal  ;
0390 }
0391
0392 //______________________________________
0393 double hcal::SignalShapeNegative(double *x, double *par)
0394 {
0395  double pedestal = par[4] + x[0]*par[5];
0396  if( x[0]<par[1]) return pedestal;
0397  double  signal = par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0398  return pedestal-signal  ;
0399 }
0400
0401
0402 //________________________________________
0403 void hcal::Save(char *filename)
0404 {
0405   fout = TFile::Open( filename, "RECREATE" );
0406   for(int ich=0; ich<NCH; ich++)
0407   {
0408     h_PulsePeaks[ich]->Write();
0409     h_tres[ich]->Write();
0410     h_ped_res[ich]->Write();
0411     h2_fit_shape[ich]->Write();
0412     //h_hilo[ich/2]->Write();
0413     //if(ich%2) h_gain[(int)ich/2]->Write();
0414     h_PulseInt[ich]->Write();
0415     h2_Pulse_H_Int[ich]->Write();
0416   }
0417   h_saturation_rates->Write();
0418   rate->Write();
0419   h_info->Write();
0420   fout->Close();
0421 }
0422