Prototype2/hcalLab/evDisplay.C

0001 #include <iomanip>
0002 #include <algorithm>
0003 #include <numeric>
0004 #include <cmath>
0005 #include <ctime>
0006
0007 #include <iostream>
0008 #include <fstream>
0009 #include <cstdlib>
0010 #include <stdint.h>
0011
0012 #include <fileEventiterator.h>
0013 #include <Event.h>
0014 #include <packet.h>
0015 #include <packet_hbd_fpgashort.h>
0016 #include <TROOT.h>
0017 #include <TF1.h>
0018 #include <TH1.h>
0019 #include <TH1F.h>
0020 #include <TCanvas.h>
0021 #include <TGraph.h>
0022 #include <TString.h>
0023 #include <TApplication.h>
0024 #include <TStyle.h>
0025 #include <TSystem.h>
0026 #include <TGClient.h>
0027 #include <TGWindow.h>
0028
0029 using namespace std;
0030
0031 //TROOT root("xxx","if you want to make a root app you can");
0032
0033 //  -----------------------------------------------------------------------------------------------
0034
0035 //const    int NCH      = 2*8;
0036 const    int NCH      = 8;        //  HG channels only
0037 //const    int NCH       = 12;
0038 const    int NSAMPLES = 24;
0039 const    int RISETIME    = 4;
0040 const    int FALLTIME    = 5;
0041 //const    int feech[]  = { 99, 98, 97, 96, 106, 107, 108, 109, 112, 113, 122, 123, 124, 125 };
0042 //const    int feech[]  = { 115, 114, 113, 112, 119, 118, 117, 116, 123, 122, 121, 120, 127, 126, 125, 124 };
0043 //const    int feech[]  = { 115, 113, 119, 117, 123, 121, 127, 125 };  //  HG channels only//  ORDERING : RIGHT/LEFT ARE TWO ENDS OF THE SAME FIBER
0044 const    int feech[]  = { 115, 119, 123, 127, 113, 117, 121, 125 };  //  HG channels only
0045
0046 TGraph  * gpulse[NCH];
0047 Double_t  shapeD[NCH][NSAMPLES];
0048 TF1     * shapes[NCH];
0049 TCanvas * ev_display, *fitFunc;
0050 TH1     * fittime, * fitmax, * fitwidth, * fitintegral, * pedval, * pedslope;
0051 TH1    ** fitPar[6];
0052 int      nevents(0);         // event counter
0053 int      channel(0);
0054
0055
0056 //  -----------------------------------------------------------------------------------------------
0057
0058 void usage()
0059 {
0060   cout << "evDisplay <prdf>" << endl;
0061 }
0062 //  -----------------------------------------------------------------------------------------------
0063 //  Unipolar pulse (0.3164*pow(x,4)*exp(-x*1.5) - Unity integral, peak ~1/3 of integral)
0064 Double_t PEDESTAL     = 2048;
0065 Double_t par0[]       = {    0., 0.,                    4., 0.,   0.,   0.};
0066 Double_t par0Max[]    = { 1000., NSAMPLES-FALLTIME+1.,  5., 1.,   2150., 0.2};
0067 Double_t par0Min[]    = {    1., RISETIME-1,            3., 2.,   1950.,-0.2};
0068 Double_t signal(0.), pedestal(0.);
0069 //  -----------------------------------------------------------------------------------------------
0070 Double_t signalShape(Double_t *x, Double_t * par){
0071   //  if(x[0]<par[1])  return PEDESTAL;
0072   pedestal = par[4]+x[0]*par[5];
0073   if(x[0]<par[1])   return pedestal;
0074   //  signal contribution
0075   signal   = par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0076   //  cout<<"EVENT  "<<nevents<<"  CHANNEL  "<<channel<<"  X  "<<x[0]<<"  signal "<<signal<<"  pedestal  "<<pedestal<<endl;
0077   return   signal+pedestal;
0078 }
0079 //  -----------------------------------------------------------------------------------------------
0080
0081 void init(){
0082   fittime     = new TH1F("fittime",    "Fitted pulse peak time",   100, 0.,   24.);
0083   fitmax      = new TH1F("fitmax",     "Fitted pulse amplitude",   500, 0.,   5000.);
0084   fitwidth    = new TH1F("fitwidth",   "Fitted pulse width",       100, 0.,   5.);
0085   pedval      = new TH1F("pedval",     "Pedestal at t==0",         100, 2000.,2100.);
0086   pedslope    = new TH1F("pedslope",   "Pedestal slope",           100, -1.,  1.);
0087   fitintegral = new TH1F("fitintegral","Integral of fitted curve (pedestal suppressed)", 100, 0., 10000.);
0088   //  fit parameters (for templating)
0089
0090
0091 }
0092 //  -----------------------------------------------------------------------------------------------
0093
0094 void fitShape(int chan){
0095   channel  = chan;
0096   TString fitName = feech[chan]%2? "HG_" : "LG_";
0097   fitName += chan;
0098   fitName += "(";
0099   fitName += feech[chan];
0100   fitName += ")";
0101   shapes[chan]= (TF1*) (gROOT->FindObject("fitName"));
0102   //  zero approximation to fit parameters
0103   //  use pulse data to find location of maximum
0104   if(shapes[chan]) delete shapes[chan];
0105   shapes[chan] = new TF1(fitName, &signalShape, 0., 24., 6);
0106   shapes[chan]->SetLineColor(2);
0107   int      peakPos = 0;
0108   Double_t peakVal = 0.;
0109   for (int iSample=0; iSample<NSAMPLES; iSample++) {
0110     if(shapeD[chan][iSample]>peakVal) {
0111       peakVal = shapeD[chan][iSample];
0112       peakPos = iSample;
0113     }
0114   }
0115   peakVal -= PEDESTAL;
0116   if(peakVal<0.) peakVal = 0.;
0117   par0[0] = peakVal/3.;
0118   par0[1] = peakPos-RISETIME;
0119   if(par0[1]<0.) par0[1] = 0.;
0120   par0[2] = 4.;
0121   par0[3] = 1.5;
0122   par0[4] = PEDESTAL;  //   we do not do pedestal subtrastion at this time
0123   par0[5] = 0;      //   slope of the pedestals is initially set to "0"
0124   shapes[chan]->SetParameters(par0);
0125   for(int parn=0; parn<6; parn++) shapes[chan]->SetParLimits(parn, par0Min[parn], par0Max[parn]);
0126   //  fitFunc->cd(chan+1);
0127   //  shapes[chan]->Draw();
0128   //  gPad->Modified();
0129   gpulse[chan]->Fit(shapes[chan], "", "RQ", 0., (Double_t)NSAMPLES);
0130   Double_t fMax = shapes[chan]->GetMaximum(par0Min[1], par0Max[1]);
0131   Double_t fPos = shapes[chan]->GetMaximumX(par0Min[1], par0Max[1]);
0132
0133
0134   cout<<"EVENT  "<<nevents<<" SHAPE  "<<(char*)(shapes[chan]->GetName())<<"  peakVal  "<<peakVal<<"  peakPos  "<<peakPos<<"  fitMax  "<<fMax<<"  fitPos  "<<fPos<<endl;
0135   ev_display->cd(chan+1);
0136   shapes[channel]->Draw("same");
0137   //  gPad->Modified();
0138
0139
0140 }
0141
0142 //  -----------------------------------------------------------------------------------------------
0143
0144 int  main(int argc, char **argv)
0145 {
0146   if ( argc!=2 )
0147     {
0148       usage();
0149       exit(-1);
0150     }
0151
0152   TString prdfname = argv[1];   // get prdf file name
0153   TString basename = prdfname;
0154   basename.ReplaceAll(".prdf","");
0155   Ssiz_t lastslash = basename.Last('/');
0156   basename.Remove(0,lastslash+1);
0157   //cout << basename << endl;
0158   basename.ReplaceAll("rc-","");
0159   //cout << basename << endl;
0160
0161   char *file = argv[1]; // get prdf file name
0162   TApplication theApp("theApp",&argc,argv);
0163
0164   //  gROOT->SetStyle("Plain");
0165 /*
0166   gStyle->SetLabelSize(0.1,"xyz");
0167   gStyle->SetTitleSize(0.08,"xyz");
0168   gStyle->SetStatH(0.4);
0169 */
0170
0171   //TCanvas *ev_display = new TCanvas("sphenix","sphenix display",800,800);
0172   //  This version of plotting is for High/Low gain preamplifier sending its outputs into two sequential readout channels
0173   int nx_c = 2;
0174   int ny_c = NCH/2+NCH%2;
0175   ev_display = new TCanvas("sphenix","sphenix display",400*nx_c,200*ny_c);
0176
0177   //  fitFunc   = new TCanvas("fitFunc","Fit Functions",  400*nx_c,200*ny_c);
0178   ev_display->SetEditable(kTRUE);
0179   //  ev_display->Divide(nx_c,ny_c);
0180   //  ev_display->SetBit(TPad::kClearAfterCR);
0181   //  fitFunc->SetEditable(kTRUE);
0182   //  fitFunc->Divide(nx_c,ny_c);
0183
0184   TString name, title;
0185   for (int ich=0; ich<NCH; ich++)
0186   {
0187     gpulse[ich] = new TGraph(NSAMPLES);
0188     name = "gch"; name += ich;
0189     gpulse[ich]->SetName(name);
0190     gpulse[ich]->SetMarkerStyle(20);
0191     gpulse[ich]->SetMarkerSize(0.4);
0192     //    gpulse[ich]->SetMarkerColor(ich+1);
0193     //    gpulse[ich]->SetLineColor(ich+1);
0194   }
0195
0196   // Open up PRDFF
0197   int status;
0198
0199   Eventiterator *it =  new fileEventiterator(file, status);
0200   if (status)
0201   {
0202     cout << "Couldn't open input file " << file << endl;
0203     delete(it);
0204     exit(1);
0205   }
0206
0207   Event *evt;
0208
0209   // Loop over events
0210   while ( (evt = it->getNextEvent()) != 0 )
0211     {
0212       int eventseq = evt->getEvtSequence();
0213
0214       if ( evt->getEvtType() != 1 ) {
0215     cout << "eventseq " << eventseq << " event type = " << evt->getEvtType() << endl;
0216     continue;
0217       }
0218       ev_display->Divide(nx_c,ny_c);
0219       ev_display->SetBit(TPad::kClearAfterCR);
0220
0221       // Get sPHENIX Packet
0222       Packet_hbd_fpgashort *spacket = dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101) );
0223       if ( spacket )
0224     {
0225       spacket->setNumSamples( NSAMPLES );
0226       int nmod_per_fem = spacket->iValue(0,"NRMODULES");
0227       //cout << "nmod_per_fem " << nmod_per_fem << endl;
0228       cout <<"EVENT "<<eventseq << "    modules per FEM "<<nmod_per_fem<<endl;
0229       nevents++;
0230       for (int ich=0; ich<NCH; ich++)
0231         {
0232           channel     = ich;
0233           for (int isamp=0; isamp<NSAMPLES; isamp++)
0234         {
0235           Short_t val = spacket->iValue(feech[ich],isamp);
0236           gpulse[ich]->SetPoint(isamp,(Double_t)isamp,(Double_t)val);
0237           shapeD[ich][isamp] = (Double_t)val;
0238         }
0239           ev_display->cd(ich+1);
0240           //gpulse[ich]->SetMaximum(4095);
0241           //gpulse[ich]->SetMinimum(0);
0242           //          if(nevents==1) gpulse[ich]->Draw("acp");
0243           gpulse[ich]->Draw("acp");
0244           //          gPad       ->Modified();
0245           fitShape(ich);
0246           //          gPad       ->Modified();
0247           //          shapes[ich]->Draw("same");
0248           //          gPad       ->Update();
0249         }
0250 //    for (int ich=0; ich<NCH; ich++)
0251 //      {
0252 //        channel    = ich;
0253 //        fitFunc    ->cd(ich+1);
0254 //        shapes[ich]->Draw();
0255 //      }
0256 //    gPad       ->Modified();
0257
0258
0259       cout<<"Updating canvas  "<<endl;
0260
0261       ev_display->Update();
0262       //      fitFunc  ->Update();
0263       //      name = basename; name += "_"; name += eventseq; name += ".png";
0264       //      ev_display->SaveAs(name);   for(int loop=0; loop<1e+6; loop++) continue;
0265       //      TGWindow * mw = (TGWindow *)(gClient->GetRoot());
0266       //      gClient->WaitFor(mw);
0267 //    if ( eventseq > 10 )
0268 //      {
0269           // string junk;
0270           // cout << "? ";
0271           // cin >> junk;
0272           // if(junk=="q") goto ALLDONE;
0273 //      }
0274
0275       delete spacket;
0276       ev_display->Clear();
0277     }
0278
0279       delete evt;
0280       //if (nevents++==10) break;
0281     }
0282
0283
0284  ALLDONE:
0285   cout<<"ALLDONE"<<endl;
0286   delete it;
0287
0288   //  theApp.Run();
0289   exit(0);
0290
0291 }
0292