Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hLabHelper.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hLabHelper.C
1 hLabHelper:: hLabHelper(){
2  // whenever it is instanciated manually (not from detector helpers) it sets all dimensions to
3  // defaults one can change to get access to something in the data
4  // if request for Helper is from one of the project specific helpers then hLabHelper job is
5  // only to handle input/output related functions. The problem - all my detector specific
6  // helpers want to create detector with all presentation tools (and there are many) first
7  // and would like to have all that material available in the root file. Unfortunately it
8  // can't be created unique before the run number is known (it's a part of its name).
9  // I want all my presentation material to be in the TFile directory so they will go outside
10  // by default and so I create a temporary storage TFile with "host" dependent name for
11  // the period the program runs and use system tools to rename that file when program is
12  // ready for termination.
13  runnumber = 0;
14  runKind = ""; // default
15  eventseq = 0;
16  eventsread = 0;
17  fhcl = NULL;
18  thcl = NULL;
19  rootTempFileName = "";
20  host = gSystem->HostName();
21  host.Remove(host.First('.'));
22 
23  if(ACTIVECHANNELS==0){
24  CHTOTAL = chinspected;
25  ACTIVECHANNELS = CHTOTAL;
26  detchannels = CHTOTAL - triggerinsp;
27  // compute total number of HG & LG channels in detector setup
28  hgDetChannels = 0;
29  lgDetChannels = 0;
30  for(int ich = 0; ich<detchannels; ich++){
31  if(feechinsp[ich]%2==HIGH) hgDetChannels ++; else lgDetChannels ++;
32  }
33  cout<<"<hLabHelper:: hLabHelper > Initializing as primary with "<<ACTIVECHANNELS<<" active channels"<<endl;
34  } else {
35  // detectorHelper will taks care of all its data - it needs only file handling utilities from labHelper
36  // create and open buffer TFile so detector helper will be able to use it to place all its presentation
37  // data in its directory
38  cout<<"<hLabHelper:: hLabHelper> Initializing as secondary";
39  rootTempFileName = (host.Contains("rcas"))? RCFrootDir : HLABrootDir;
40  TDatime T;
41  TString fN = "hcalHelper_"; fN += T.GetDate(); fN += "_"; fN += T.GetTime(); fN += ".root";
42  rootTempFileName += fN;
43  cout<<"<hLabHelper:: hLabHelper> ROOT graphic objects will be written into TFile "<<rootTempFileName<<endl;
44  fhcl = new TFile(rootTempFileName, "recreate");
45  FileStat_t buf;
46  if(gSystem->GetPathInfo(rootTempFileName, buf)) {
47  cout<<"<runToRootFile> : File "<<rootTempFileName<<" still does not exist"<<endl;
48  } else {cout<<"<hLabHelper:: hLabHelper> ROOT file "<<rootTempFileName<<" created"<<endl;}
49  fhcl->cd();
50  return;
51  }
52  active = new Int_t [ACTIVECHANNELS];
53  adc = new Float_t * [ACTIVECHANNELS];
54  fitResults = new Double_t * [ACTIVECHANNELS];
55  for (int ich=0; ich<ACTIVECHANNELS; ich++) {
56  adc[ich] = new Float_t [NSAMPLES];
57  fitResults[ich] = new Double_t [NPARAMETERS];
58  }
59  pedestal = new Float_t [ACTIVECHANNELS]; // [ACTIVECHANNELS] */
60  rawPeak = new Float_t [ACTIVECHANNELS];
61  rawTime = new Float_t [ACTIVECHANNELS];
62  fitPeak = new Float_t [ACTIVECHANNELS];
63  fitTime = new Float_t [ACTIVECHANNELS];
64  fitInt = new Float_t [ACTIVECHANNELS];
65  fitChi2 = new Float_t [ACTIVECHANNELS];
66  fitPeakRMS = new Float_t [ACTIVECHANNELS];
67  fitTimeRMS = new Float_t [ACTIVECHANNELS];
68  rdata = new TH2F("rdata", "Raw Peak Values for all active channels",ACTIVECHANNELS, 0, ACTIVECHANNELS, 1024, -2048., 2048.);
69  shapes = new TF1 * [ACTIVECHANNELS];
70  sigFit = new TF1 * [ACTIVECHANNELS];
71  rvsc = new TF1 * [ACTIVECHANNELS];
72  fmsc = new TF1 * [ACTIVECHANNELS];
73  gpulse = new TGraph * [ACTIVECHANNELS];
74  ft = new TH1 * [ACTIVECHANNELS];
75  fm = new TH1 * [ACTIVECHANNELS];
76  fw = new TH1 * [ACTIVECHANNELS];
77  fint = new TH1 * [ACTIVECHANNELS];
78  fpd = new TH1 * [ACTIVECHANNELS];
79  fchi2 = new TH1 * [ACTIVECHANNELS];
80  fitPar = new TH1 ** [ACTIVECHANNELS];
81  rpeak = new TH1 * [ACTIVECHANNELS];
82  rtm = new TH1 * [ACTIVECHANNELS];
83  for (int ich = 0; ich < ACTIVECHANNELS; ich++) {
84  fitPar[ich] = new TH1 * [NPARAMETERS];
85  shapes[ich] = NULL;
86  gpulse[ich] = NULL;
87  }
88  for (int ig = 0; ig<2; ig++){
89  evtadcsum[ig] = new Double_t [NSAMPLES];
90  evtGraph[ig] = new TGraph;
91  TString sum = "adcsum_"; sum += (ig==HIGH)? "HG" : "LG";
92  TString sig = "signal_"; sig += (ig==HIGH)? "HG" : "LG";
93  evtShape[ig] = new TF1(sum, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
94  evtSignal[ig] = new TF1(sig, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
95  // Double_t times[ACTIVECHANNELS];
96  // Double_t chi2 [ACTIVECHANNELS];
97  evtfitpar[ig] = new Double_t [NPARAMETERS];
98  }
99  }
100 // **************************************************************************
101 
102 bool hLabHelper::evDisplay(Int_t rn){
103  if(rn<=0) return false;
104  runnumber = rn;
105  if(!runToPRDFFile(rn)) {
106  cout<<"<setPRDFRun> Failed to find the .prdf file for run number "<<rn<<endl;
107  return false;
108  }
109  // Open up PRDFF
110  int status;
111  Eventiterator *it = new fileEventiterator(prdfName, status);
112  if (status)
113  {
114  cout << "Couldn't open input file " << prdfName << endl;
115  delete(it);
116  return false;
117  }
118  // update channel map
119  updateMap();
120  int chperpanel = DISPLAYX*DISPLAYX;
121  int panels = (ACTIVECHANNELS%chperpanel==0? ACTIVECHANNELS/chperpanel : ACTIVECHANNELS/chperpanel+1);
122  int nx_c = DISPLAYX;
123  int ny_c(0);
124  if(panels>1||ACTIVECHANNELS%chperpanel==0) {
125  ny_c = DISPLAYX;
126  } else {
127  ny_c = ACTIVECHANNELS/DISPLAYX + (ACTIVECHANNELS%DISPLAYX==0? 0 : 1);
128  }
129 
130  // cout<<"Display "<<panels<<" "<<nx_c<<"/"<<ny_c<<endl;
131 
132  gStyle->SetOptFit(111);
133  gStyle->SetFitFormat("5.2g");
134  // creating Canvases
135  TCanvas * evdisplay[panels];
136  TString cnvsName = "EvDisplay_";
137  TString cnvsTitle = "sPhenix Channel display for CHANNELS ";
138  for (int iC = 0; iC < panels; iC++) {
139  TString cN = cnvsName;
140  cN += iC;
141  TString cT = cnvsTitle;
142  cT += (chperpanel*iC);
143  cT += " - ";
144  cT += (chperpanel*(iC+1)-1);
145  evdisplay[iC] = new TCanvas(cN, cT, 300*nx_c, 150*ny_c);
146  evdisplay[iC]->SetEditable(kTRUE);
147  evdisplay[iC]->Divide(nx_c,ny_c);
148  // if(iC==panels) {
149  // evdisplay[iC]->SetBit(TPad::kClearAfterCR);
150  // }
151  }
152  // ====================================================================================
153  // now TGraph for every channel
154  for (int ich = 0; ich<ACTIVECHANNELS; ich++)
155  {
156  gpulse[ich] = new TGraph(NSAMPLES);
157  TString name = "gch"; name += ich;
158  gpulse[ich]->SetName(name);
159  name = "Channel "; name += ich; name += " / "; name += feechinsp[ich]; name += (feechinsp[ich]%2==HIGH? " HG" : " LG");
160  gpulse[ich]->SetTitle(name);
161  gpulse[ich]->SetMarkerStyle(20);
162  gpulse[ich]->SetMarkerSize(0.4);
163  //gpulse[ich]->SetMarkerColor(ich+1);
164  //gpulse[ich]->SetLineColor(ich+1);
165  }
166  Event *evt;
167 
168  // Loop over events
169  while ( (evt = it->getNextEvent()) != 0 )
170  {
171  eventseq = evt->getEvtSequence();
172 
173  if ( evt->getEvtType() != 1 ) {
174  cout << "eventseq " << eventseq << " event type = " << evt->getEvtType() <<" ignored "<< endl;
175  } else {
176  // Get sPHENIX Packet
177  Packet_hbd_fpgashort *spacket = dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101) );
178  if ( spacket ) {
179  spacket->setNumSamples( NSAMPLES );
180  cout <<"PRDF FILE "<<prdfName<<" EVENT "<<eventseq <<endl;
181  for (int ich=0; ich<ACTIVECHANNELS; ich++) {
182  for (int is=0; is<NSAMPLES; is++) {
183  adc[ich][is] = spacket->iValue(active[ich],is);
184  // cout<<"<evDisplay Event "<< eventseq<<" ich/is/adc "<<ich<<" - "<<is<<" - "<<adc[ich][is]<<endl;
185  }
186  }
187  dofixes();
188  getDetectorTiming();
189  // cout<<"<TileFit> Peak "<< evtsum.evtpeak<<" Time "<<evtsum.evttime<<" Chi2 "<<evtsum.tChi2 <<endl;
190  displaySumEvent();
191  // we do generic display of square TCanvases DISPLAYX x DISPLAYX size
192  // for (int evD = 0; evD<panels; evD++){
193  // TString evDName = cnvsName + evD;
194  // TString evDTitle = cnvsTitle + evD;
195  // evdisplay->SetBit(TPad::kClearAfterCR);
196  // evdisplay->Divide(nx_c,ny_c);
197 
198  for (int ich=0; ich<ACTIVECHANNELS; ich++) {
199  for (int is=0; is<NSAMPLES; is++) gpulse[ich]->SetPoint(is,(Double_t)is, adc[ich][is]);
200  int iC = ich/chperpanel;
201  int iPad = (ich-chperpanel*iC);
202  int iPx = iPad%nx_c;
203  int iPy = ny_c - iPad/nx_c - 1;
204  iPad = iPy*nx_c + iPx + 1;
205  evdisplay[iC] ->cd(iPad);
206  // cout<<iC<< " "<<ich<<" "<<iPx<<" "<<iPy<<" "<<iPad<<endl;
207  gpulse[ich]->Draw("acp");
208  // Double_t peakVal;
209  // Int_t peakPos;
210  // fitShape(ich, peakVal, peakPos, (ich<detchannels? 1 : 0));
211  // Double_t fMax = shapes[ich]->GetMaximum(par0Min[1], par0Max[1]);
212  // Double_t fPos = shapes[ich]->GetMaximumX(par0Min[1], par0Max[1]);
213  // Double_t fChi2 = shapes[ich]->GetChisquare()/shapes[ich]->GetNDF();
214  // cout<<"EVENT "<<eventseq<<" SHAPE "<<(char*)(shapes[ich]->GetName())<<" peakVal "<<peakVal<<" peakPos "<<peakPos<<" fitMax "<<fMax<<" fitPos "<<fPos<<" fChi2 "<<fChi2<<endl;
215  // shapes[ich]->Draw("same");
216  // evdisplay[iC]->Update();
217  }
218  for (int iC = 0; iC < panels; iC++) evdisplay[iC]->Update();
219  char dummy;
220  cin>>dummy;
221  // }
222  delete spacket;
223  // for (int iC = 0; iC < panels; iC++) evdisplay[iC]->Clear();
224  } else cout<<"<evDisplay> Packet hbd not found. EventSeq "<<eventseq<<" EventType "<< evt->getEvtType()<<endl;
225  }
226  delete evt;
227  }
228  cout<<"<evDisplay> EoF reached "<<endl;
229 
230  return true;
231 }
232 
233 // **************************************************************************
234 void hLabHelper::updateMap(){
235  for (int ich = 0; ich<ACTIVECHANNELS; ich++) active[ich] = feechinsp[ich] ;
236 }
237 
238 // **************************************************************************
239 void hLabHelper::setRunKind(TString kind){
240  runKind = kind;
241  cout<<"hLabHelper::setRunKind> will look for run data in "<<runKind<<endl;
242 }
243 
244 
245 
246 // **************************************************************************
247 void hLabHelper::getFileLists(){
248  // cout<<"<getFileLists> hostname "<<host<<endl;
249  // ====================================================================================
250  // parsing .root file name from runnumber
251  rootdirname = (host.Contains("rcas"))? RCFrootDir : HLABrootDir;
252  prdfdirname = (host.Contains("rcas"))? RCFdataDir : HLABdataDir;
253  // in t1044 we have "junk", "cosmics", "led" and "beam" kinds of events
254  cout<<"<hLabHelper::getFileLists> "<<rootdirname<<endl;
255  cout<<"<hLabHelper::getFileLists> "<<prdfdirname<<endl;
256  if(!runKind.IsNull()) {
257  prdfdirname += runKind; prdfdirname +="/";
258  cout<<"hLabHelper::getFileLists> PRDF Directory name updated "<<prdfdirname<<endl;
259  }
260  TSystemDirectory rdir(rootdirname, rootdirname);
261  TSystemDirectory pdir(prdfdirname, prdfdirname);
262  TList * rootfiles = rdir.GetListOfFiles();
263  TList * prdffiles = pdir.GetListOfFiles();
264  Int_t rn;
265  TString rid;
266  if(rootfiles){
267  TSystemFile * file;
268  TString fname;
269  TIter next(rootfiles);
270  while((file = (TSystemFile *) next())) {
271  fname = file->GetName();
272  if(!file->IsDirectory() && fname.EndsWith(".root")) {
273  // .root file found - get run number
274  decoderun(fname, rn, rid);
275  TString name = rootdirname;
276  name += fname;
277  rootfile rf(name, rn, rid);
278  roots.push_back(rf);
279  // cout<<".root entry "<<name<<endl;
280  } else {
281  // cout<<"<getFileList> candidate root "<<fname<<" ignored "<<endl;
282 
283  }
284  }
285  }
286  if(prdffiles){
287  TSystemFile * file;
288  TString fname;
289  TIter next(prdffiles);
290  while((file = (TSystemFile *) next())) {
291  fname = file->GetName();
292  // cout<<fname<<endl;
293  if(!file->IsDirectory() && fname.EndsWith(".prdf")) {
294  // .prdf file found - get run number
295  decoderun(fname, rn, rid);
296  TString name = prdfdirname;
297  name += fname;
298  prdffile rf(name, rn, rid);
299  prdfs.push_back(rf);
300  // cout<<".prdf entry "<<name<<endl;
301  } else {
302  // cout<<"<getFileList> candidate prdf "<<fname<<" ignored "<<endl;
303 
304  }
305  }
306 
307  }
308  cout<<"<getFileLists> .root files found "<<roots.size()<<" .prdf files found "<<prdfs.size()<<endl;
309 }
310 
311 // **************************************************************************
312 void hLabHelper::decoderun(TString & fname, Int_t & rn, TString & rid){
313  if(fname.Contains(".root") || fname.Contains(".prdf")) {
314  // cout<<"<decoderun> : HCAL file "<<fname<<" confirmed"<<endl;
315  rid = fname; // get run number from root file name
316  rid.ReplaceAll(".root",""); rid.ReplaceAll(".prdf","");
317  Ssiz_t lastslash = rid.Last('/');
318  rid.Remove(0,lastslash+1);
319  rid.ReplaceAll("HClab_",""); rid.ReplaceAll("rc-",""); rid.ReplaceAll("cosmics_",""); rid.ReplaceAll("beam_","");
320 
321  TString id = rid;
322  id.Remove(rid.Last('-'));
323  rn = id.Atoi();
324  // cout<<"<decoderun> "<<fname<<" runnumber "<<rn<<" runId "<<rid<<endl;
325  }
326 }
327 
328 // **************************************************************************
329 TFile * hLabHelper::attachrootrun(Int_t rn){
330  if(!roots.size()) {
331  cout<<"<attachrootrun> No .root files in the list: Collecting "<<endl;
332  getFileLists();
333  }
334  // find the rootfile rn in the local list of files kept by Helper
335  std::list<rootfile>::iterator it = roots.begin();
336 
337  while (it!=roots.end()) {
338  // it->print();
339  if(it->runnumber==rn) {
340  runnumber = rn;
341  runId = it->runid;
342  cout<<"<hLabHelper::attachrootrun> ROOT file with "<<runId<<" identifier found: "<<it->name<<endl;
343  // this is dangerous
344  return new TFile(it->name);
345  } else it++;
346  }
347  cout<<"<hLabHelper::attachrootrun> Requested ROOT file for Run "<<rn<<" not found"<<endl;
348  return NULL;
349  // auto it = std::find_if( std::begin( roots ),
350  // std::end( roots ),
351  // [&]( const rootfile rf ){ return 1 == ( rf.runnumber==rn ); } );
352  // if(roots.end() == it) return NULL;
353  // return new TFile(it->name);
354 }
355 
356 
357 
358 // const int pos = std::distance( myList.begin(), it ) + 1;
359 // std::cout << "item is found at position " << pos << std::endl;
360 
361 
362 // **************************************************************************
363 Bool_t hLabHelper::setPRDFRun(int run, Bool_t bookH){
364  if(run<=0) return false;
365  runnumber = run;
366  if(!runToPRDFFile(run)) {
367  cout<<"<setPRDFRun> Failed to find the .prdf file for run number "<<run<<endl;
368  return false;
369  }
370  initPRDFRun(bookH);
371  return true;
372 }
373 
374 
375 
376 // **************************************************************************
377 Bool_t hLabHelper::setROOTRun(int run){
378  if(!runToRootFile(run)) {
379  cout<<"<setROOTRun> Failed to find the .root file for run number "<<run<<endl;
380  return false;
381  }
382  return true;
383 }
384 
385 
386 // **************************************************************************
387 
388 // The assumption hear is that the .root file is either already attached to this ROOT session or
389 // must be looked for using its runnumber
390 Int_t hLabHelper:: runToRootFile(Int_t run){
391  if(run<=0) {
392  cout<<"<runToRootFile> WARNING : Requested run number = "<<run<<endl;
393  TSeqCollection * sq = gROOT->GetListOfFiles();
394  if(!sq->IsEmpty()) {
395  // there are already attached files - let's see if we may find something over there
396  // now it is really .rootf into runId not what is in the name
397  TString current = ((TObject *)(sq->First()))->GetName();
398  // check the extension first
399  if(current.Contains(".root")) {
400  cout<<"<hLabHelper::runToRootFile> : attached .root file "<<current<<" found"<<endl;
401  rootfName = current;
402  runId = current; // get run number from root file name
403  runId.ReplaceAll(".root","");
404  Ssiz_t lastslash = runId.Last('/');
405  runId.Remove(0,lastslash+1);
406  runId.ReplaceAll("HClab_","");
407  TString rn = runId;
408  rn.Remove(rn.Last('-'));
409  runnumber = rn.Atoi();
410  return(runnumber);
411  }
412  } else return 0;
413  } else {
414  // TODO: check if TFile for requested run is already attached
415  // let's try to find and attach file with the right runnumber
416  // find where I am
417  cout<<"<runToRootFile> : looking into .root directory for runnumber "<<run<<endl;
418  // ====================================================================================
419  // parsing .root file name from runnumber
420  rootfName = (host.Contains("rcas"))? RCFrootDir : HLABrootDir;
421  TString fname = rootfName;
422  FileStat_t buf;
423  runId = "00"; runId += run; runId += "-0";
424  fname += "HClab_"; fname += runId; fname += ".root";
425  if(gSystem->GetPathInfo(fname, buf)) {
426  fname = rootfName;
427  runId = "0"; runId += run; runId += "-0";
428  fname += "HClab_"; fname += runId; fname += ".root";
429  rootfName = fname;
430  if(gSystem->GetPathInfo(fname, buf)) {
431  cout<<"<runToRootFile> : File "<<fname<<" does not exist"<<endl;
432  } else cout<<"<runToRootFile> : File "<<fname<<" already exist and will be recreated"<<endl;
433  }
434  runnumber = run;
435  }
436  if(!rootTempFileName.IsNull()) {
437  cout<<"<hLabHelper:: runToRootFile> found temporaly root file "<<rootTempFileName<<endl; } else {
438  fhcl = new TFile(rootfName,"recreate");
439  cout<<"hLabHelper:: runToRootFile> "<<rootfName<<" opened"<<endl;
440  }
441  return runnumber;
442 }
443 
444 // **************************************************************************
445 
446 Int_t hLabHelper::runToPRDFFile(Int_t run){
447  // ====================================================================================
448  // find where I am
449  // ====================================================================================
450  // parsing .prdf file name from runnumber
451  prdfName = (host.Contains("rcas"))? RCFdataDir : HLABdataDir;
452  if(!runKind.IsNull()) {
453  prdfName += runKind; prdfName += "/";
454  cout<<"hLabHelper::runToPRDFFile> PRDF Directory name updated "<<prdfName<<endl;
455  }
456  TString fname = prdfName;
457  FileStat_t buf;
458  // runId = "00"; runId += run; runId += "-0";
459  // fname += "rc-"; fname += runId; fname += ".prdf";
460  runId = "0000"; runId += run; runId += "-0000";
461  fname += runKind; fname += "_";
462  fname += runId; fname += ".prdf";
463  if(gSystem->GetPathInfo(fname, buf)) {
464  fname = prdfName;
465  // runId = "00"; runId += run; runId += "-0";
466  // fname += "rc-"; fname += runId; fname += ".prdf";
467  runId = "0000"; runId += runnumber; runId += "-0000";
468  fname += runKind; fname += "_";
469  fname += runId; fname += ".prdf";
470  if(gSystem->GetPathInfo(fname, buf)) {
471  cout<<"<runToPRDFFile> File "<<fname<<" with runId "<<runId<<" does not exist"<<endl;
472  return 0;
473  }
474  }
475  prdfName = fname;
476  // ====================================================================================
477  cout<<"<runToPRDFFile> Parsed prdf file name "<<prdfName<<endl;
478  // ====================================================================================
479 
480  cout << "<runToPRDFFile> prdfName "<<prdfName << " RunId "<<runId<<" Run Number "<<runnumber<<endl;
481 
482  // ====================================================================================
483  return runnumber;
484 }
485 
486 
487 // **************************************************************************
488 void hLabHelper::initPRDFRun(Bool_t bookH){
489  runToRootFile(runnumber);
490  // ====================================================================================
491  // initialize root file and TTree (the idea is to save all data and histos at the end of run)
492  if(runnumber<=0) {
493  cout<<"<hLabHelper::initPRDFRun> Cant instantiate run with runnumber = 0. Exiting"<<endl;
494  return;
495  }
496  // ====================================================================================
497  // TApplication theApp("theApp",&argc,argv);
498 
499  gROOT->SetStyle("Plain");
500  /*
501  gStyle->SetLabelSize(0.1,"xyz");
502  gStyle->SetTitleSize(0.08,"xyz");
503  gStyle->SetStatH(0.4);
504  */
505  if(!bookH) return;
506  // fit parameters (for templating)
507  for (int ich = 0; ich<ACTIVECHANNELS; ich ++){
508  TString fhn = "ft_"; fhn += ich; TString fht = "Fitted pulse peak time CH="; fht += ich;
509  ft[ich] = new TH1F(fhn, fht, 100, 0., NSAMPLES);
510  fhn = "fm_"; fhn += ich; fht = "Fitted pulse amplitude CH="; fht += ich;
511  fm[ich] = new TH1F(fhn, fht, 4100, 0., 4100.);
512  fhn = "fw_"; fhn += ich; fht = "Fitted pulse width CH="; fht += ich;
513  fw[ich] = new TH1F(fhn, fht, 100, 0., 5.);
514  fhn = "fpd_"; fhn += ich; fht = "Fitted pedestal at t=PEAK CH="; fht += ich;
515  fpd[ich] = new TH1F(fhn, fht, 200, 2000., 2100.);
516  fhn = "fint_"; fhn += ich; fht = "Integral of fitted curve (pedestal suppressed) CH="; fht += ich;
517  fint[ich] = new TH1F(fhn, fht, 4000, 0., 4000.);
518  fhn = "fchi2_";fhn += ich; fht = "Chi2 of fit CH="; fht += ich;
519  fchi2[ich] = new TH1F(fhn, fht, 500, 0., 100.);
520 
521  for(int ip=0; ip<NPARAMETERS; ip++){ TString phn = "p_"; phn += ip; phn += "_"; phn += ich;
522  TString pht = "Fit Parameter "; pht += ip; pht += " Channel "; pht += ich;
523  fitPar[ich][ip] = new TH1F(phn, pht, 100, par0Min[ip], par0Max[ip]);
524  }
525 
526  TString rvn = "rv_"; rvn += ich; TString rvt = "Raw peak value CH = "; rvt += ich;
527  rpeak[ich] = new TH1F(rvn, rvt, 4100, 0., 4100.);
528  TString rtn = "rt_"; rtn += ich; TString rtt = "Raw peak time CH = "; rtt += ich;
529  rtm[ich] = new TH1F(rtn, rtt, NSAMPLES, 0., NSAMPLES);
530  }
531  for (int ich=0; ich<ACTIVECHANNELS; ich++) {
532  TString fn = "S_"; fn += ich;
533  sigFit[ich] = new TF1(fn, &signalShape, 0., (Double_t)NSAMPLES, 6);
534  // sigFit[ich]->AbsValue(true);
535  } // This version of plotting is for High/Low gain preamplifier sending its outputs into two sequential readout channels
536  // int nx_c = 2;
537  // int ny_c = ACTIVECHANNELS/2+ACTIVECHANNELS%2;
538  // ev_display = new TCanvas("sphenix","sphenix display",400*nx_c,200*ny_c);
539 
540  // fitFunc = new TCanvas("fitFunc","Fit Functions", 400*nx_c,200*ny_c);
541  // ev_display->SetEditable(kTRUE);
542  // ev_display->Divide(nx_c,ny_c);
543  // ev_display->SetBit(TPad::kClearAfterCR);
544  // fitFunc->SetEditable(kTRUE);
545  // fitFunc->Divide(nx_c,ny_c);
546 
547  TString name, title;
548  for (int ich=0; ich<ACTIVECHANNELS; ich++)
549  {
550  gpulse[ich] = new TGraph(NSAMPLES);
551  name = "gch"; name += ich;
552  gpulse[ich]->SetName(name);
553  gpulse[ich]->SetMarkerStyle(20);
554  gpulse[ich]->SetMarkerSize(0.4);
555  //gpulse[ich]->SetMarkerColor(ich+1);
556  //gpulse[ich]->SetLineColor(ich+1);
557  }
558 
559 }
560 
561 // -----------------------------------------------------------------------------------------------
562 void hLabHelper::renameAndCloseRF(){
563  if(rootTempFileName.IsNull()) return;
564  fhcl -> Close();
565  gSystem->Rename(rootTempFileName, rootfName);
566  cout<<"<hLabHelper::copyAndCloseRootFile> Renamed "<<endl;
567  cout<<" "<<rootTempFileName<<endl;
568  cout<<"into "<<endl;
569  cout<<" "<<rootfName<<endl;
570 }
571 
572 
573 // -----------------------------------------------------------------------------------------------
574 void hLabHelper::makeCanvasDirectory(){
575  rootdirname = (host.Contains("rcas"))? RCFrootDir : HLABrootDir;
576  cDirectory = rootdirname; cDirectory += "canvases/r";
577  cDirectory += runnumber; cDirectory += "/";
578  FileStat_t buf;
579  if(!gSystem->GetPathInfo(cDirectory, buf)) {
580  cout<<"<hLabHelper::makeCanvasDirectory> Directory "<<cDirectory<<" already exists "<<endl;
581  } else {
582  // gSystem->MakeDirectory(cDirectory);
583  TString command = "mkdir "; command += cDirectory;
584  gSystem->Exec(command);
585  cout<<"<hLabHelper::makeCanvasDirectory> Directory "<<cDirectory<<" created "<<endl;
586 
587  }
588 }
589 // -----------------------------------------------------------------------------------------------
590 
591 void hLabHelper::fitShape(int ich, Double_t & rVal, Int_t & rTime, Int_t mode){
592  if(!shapes[ich]){
593  TString fitName = active[ich]%2? "HG_" : "LG_";
594  fitName += ich/2;
595  fitName += "(";
596  fitName += ich;
597  fitName += ")";
598  shapes[ich] = new TF1(fitName, &signalShape, 0., (Double_t)NSAMPLES, 6);
599  }
600  // GAIN range
601  Int_t ig = feechinsp[ich]%2==HIGH? HIGH : LOW;
602  // zero approximation to fit parameters
603  // use pulse data to find location of maximum
604  // cout<<"<fitShape> Event "<<hclt->eventseq<<" channel "<<ich<<endl;
605  rTime = 0;
606  rVal = 0.;
607  Int_t iss(RISETIME-2), ise(NSAMPLES-FALLTIME);
608  Int_t rMaxVal = adc[ich][iss]; Int_t maxsmpl = iss;
609  Int_t rMinVal = adc[ich][iss]; Int_t minsmpl = iss;
610 
611  for (int is = iss; is<=ise; is++) {
612  if(adc[ich][is]>rMaxVal) {
613  rMaxVal = adc[ich][is];
614  maxsmpl = is;
615  }
616  if(adc[ich][is]<rMinVal) {
617  rMinVal = adc[ich][is];
618  minsmpl = is;
619  }
620  // cout<<ich<<" "<<is<<" "<<adc[ich][is]<<" "<<rVal<<" "<<rTime<<endl;
621  }
622  rMinVal -= PEDESTAL; rMaxVal -= PEDESTAL;
623  if(abs(rMinVal) > abs(rMaxVal)) {
624  rVal = rMinVal; rTime = minsmpl;
625  } else {
626  rVal = rMaxVal; rTime = maxsmpl;
627  }
628 
629  // cout<<" rVal-ped "<<rVal<<endl;
630 
631  par0[0] = rVal/3.;
632  if(mode == 0 ) {
633  par0[1] = max(0.,(Double_t)(rTime-RISETIME));
634  par0[2] = 4.;
635  par0[3] = 1.6;
636  par0[4] = PEDESTAL;
637  // we do not do pedestal subtrastion at this time
638  par0[5] = 0.; // slope of the pedestals is initially set to "0"
639  shapes[ich]->SetParameters(par0);
640  for(int parn=0; parn<6; parn++) shapes[ich]->SetParLimits(parn, par0Min[parn], par0Max[parn]);
641  } else {
642  // channels which belong to detector
643  par0[1] = evtfitpar[ig][1];
644  par0[2] = evtfitpar[ig][2];
645  par0[3] = evtfitpar[ig][3];
646  par0[4] = PEDESTAL;
647  par0[5] = 0.;
648  shapes[ich]->SetParameters(par0);
649  shapes[ich]->SetParLimits(2, evtfitpar[ig][2], evtfitpar[ig][2]);
650  shapes[ich]->SetParLimits(3, evtfitpar[ig][3], evtfitpar[ig][3]);
651  }
652  // fitFunc->cd(ich+1); // shapes[ich]->Draw();
653  // gPad->Modified();
654  gpulse[ich]->Fit(shapes[ich], "RQ", "NQ", 0., (Double_t)NSAMPLES);
655  shapes[ich]->GetParameters(fitResults[ich]);
656  // ev_display->cd(ich+1);
657  // shapes[hHelper->ich]->Draw("same");
658 
659 }
660 
661 
662 // **************************************************************************
663 // The next few functions are to process data from HCal Lab recorded via HBD electronics
664 
665 void hLabHelper::collect(Bool_t fitshape){
666  // Fit parameters
667 
668  // Double_t uSum(0.), cSum(0.), luSum(0.), lcSum(0.), ruSum(0.), rcSum(0.);
669  // Double_t YU(0.), YC(0.), XU(0.), XC(0.);
670  // hits above threshold in this event depending on the threshold (binned in Chi units - 8 counts)
671  Double_t peakVal;
672  Int_t peakPos;
673  for (int ich=0; ich<ACTIVECHANNELS; ich++) {
674  for (int is=0; is<NSAMPLES; is++){
675  gpulse[ich]->SetPoint(is,(Double_t)is, adc[ich][is]);
676  }
677  if(fitshape) fitShape(ich, peakVal, peakPos, (ich<detchannels? 1 : 0));
678  // fitShape(ich, peakVal, peakPos, 1);
679  }
680 
681  for (int ich=0; ich<ACTIVECHANNELS; ich++){
682  int rTime = 0;
683  Double_t rVal = 0.;
684  // Int_t iss(RISETIME-2), ise(NSAMPLES-FALLTIME);
685  Int_t iss(0), ise(NSAMPLES-1);
686  Int_t rMaxVal = adc[ich][iss]; Int_t maxsmpl = iss;
687  Int_t rMinVal = adc[ich][iss]; Int_t minsmpl = iss;
688 
689  for (int is = iss; is<=ise; is++) {
690  if(adc[ich][is]>rMaxVal) {
691  rMaxVal = adc[ich][is];
692  maxsmpl = is;
693  }
694  if(adc[ich][is]<rMinVal) {
695  rMinVal = adc[ich][is];
696  minsmpl = is;
697  }
698  }
699  rMinVal -= PEDESTAL; rMaxVal -= PEDESTAL;
700  if(abs(rMinVal) > abs(rMaxVal)) {
701  rVal = rMinVal; rTime = minsmpl;
702  } else {
703  rVal = rMaxVal; rTime = maxsmpl;
704  }
705  // cout<<"<collect> "<<ich<<" "<<rTime<<" "<<adc[ich][(int)rTime]<<endl;
706  rawTime[ich] = rTime;
707  rawPeak[ich] = rVal;
708  rawPeak[ich] -= (fitshape? (fitResults[ich][4]+fitResults[ich][5]*rTime) : PEDESTAL); //
709  rpeak[ich] -> Fill(rawPeak[ich]); rtm[ich] -> Fill(rawTime[ich]);
710  rdata -> Fill(ich, rawPeak[ich]);
711  // Below only in case if fitshape is required
712  if(fitshape){
713  sigFit[ich]->SetParameters(fitResults[ich]);
714  sigFit[ich]->SetParameter(4, 0.);
715  sigFit[ich]->SetParameter(5, 0.);
716 
717  // Double_t pVal = sigFit[ich]->GetMaximum(par0Min[1], par0Max[1]);
718  // Double_t pTime = sigFit[ich]->GetMaximumX(par0Min[1], par0Max[1]);
719  Double_t sgp = sigFit[ich]->GetMaximum(par0Min[1], par0Max[1]);
720  Double_t tmp = sigFit[ich]->GetMaximumX(par0Min[1], par0Max[1]);
721  Double_t sgn = sigFit[ich]->GetMinimum(par0Min[1], par0Max[1]);
722  Double_t tmn = sigFit[ich]->GetMinimumX(par0Min[1], par0Max[1]);
723  Double_t pVal = rVal;
724  Double_t pTime = rTime;
725  if(tmp!=0&&abs(sgp)>abs(sgn)) {
726  pVal = sgp;
727  pTime = tmp;
728  // cout<<"<getDetectorTimingT>: Positive peak selected"<<endl;
729  } else {
730  pVal = sgn;
731  pTime = tmn;
732  // cout<<"<getDetectorTimingT>: Negative peak selected"<<endl;
733  }
734  Double_t xmin = max(0.,pTime-RISETIME);
735  Double_t xmax = min(Double_t(NSAMPLES), pTime+FALLTIME);
736  Double_t pInt = (pVal!=0? sigFit[ich]->Integral(xmin, xmax) : 0);
737  // if(pInt==0) {
738  // cout<<"<collect> Integral zero Event "<<hclt->eventseq<<" Ch "<<ich<<" pVal/pTime/xmin/xmax "<<pVal<<" "<<pTime<<" "<<xmin<<" "<<xmax<<endl;
739  // }
740  // pVal -= fped;
741  fitPeak[ich] = pVal;
742  fitTime[ich] = pTime;
743  fitInt[ich] = pInt;
744  fitChi2[ich] = shapes[ich]->GetChisquare()/shapes[ich]->GetNDF();
745  // trigger contribution
746  ft[ich] -> Fill(pTime);
747  fm[ich] -> Fill(pVal);
748  // Double_t mean = sigFit[ich]->Mean(xmin, xmax);
749  Double_t mom2 = (pVal!=0? sigFit[ich]->CentralMoment(2, xmin, xmax) : 0.);
750  // if(mom2==0) {
751  // cout<<"<collect> CentralMoment zero Event "<<hclt->eventseq<<" Ch "<<ich<<" pVal/pTime/pInt/xmin/xmax "<<pVal<<" "<<pTime<<" "<<pInt<<" "<<xmin<<" "<<xmax<<endl;
752  // }
753  Double_t rms = sqrt(mom2);
754  fw[ich] -> Fill(rms);
755  fpd[ich] -> Fill(fitResults[ich][4]+pTime*fitResults[ich][5]);
756  fint[ich] -> Fill(pInt);
757  fchi2[ich] -> Fill(fitChi2[ich]);
758 
759  for (int ip=0; ip<NPARAMETERS; ip++) fitPar[ich][ip]->Fill(fitResults[ich][ip]);
760  // Double_t fVal = GetMaximum( par0Min[2], par0Max[2]);
761  }
762  }
763 
764  if(!reject()){
765  // few histogramms for clean sample
766  };
767 }
768 
769 // **************************************************************************
770 
771 void hLabHelper::getDetectorTiming(){
772  // build eventsum (adc sum without pedestal subtruction)
773  // do fitting etc (HG/LG separately)
774 
775  for(int is = 0; is < NSAMPLES; is++){
776  evtadcsum[0][is] = 0.;
777  evtadcsum[1][is] = 0.;
778  for (int ich = 0; ich < detchannels; ich++) {
779  Int_t ig = feechinsp[ich]%2==HIGH? HIGH : LOW;
780  evtadcsum[ig][is] += adc[ich][is];
781  }
782  }
783  for (int is = 0; is < NSAMPLES; is++) {
784  evtGraph[HIGH]->SetPoint(is,(Double_t)is,evtadcsum[HIGH][is]);
785  evtGraph[LOW] ->SetPoint(is,(Double_t)is,evtadcsum[LOW][is]);
786  }
787  Double_t rVal, rTime;
788  for(int ig = 0; ig<2; ig++) {
789  Int_t iss(RISETIME-2), ise(NSAMPLES-FALLTIME);
790  Int_t rMaxVal = evtadcsum[ig][iss]; Int_t maxsmpl = iss;
791  Int_t rMinVal = evtadcsum[ig][iss]; Int_t minsmpl = iss;
792 
793  for (int is = iss; is<=ise; is++) {
794  if(evtadcsum[ig][is]>rMaxVal) {
795  rMaxVal = evtadcsum[ig][is];
796  maxsmpl = is;
797  }
798  if(evtadcsum[ig][is]<rMinVal) {
799  rMinVal = evtadcsum[ig][is];
800  minsmpl = is;
801  }
802  // cout<<ich<<" "<<is<<" "<<adc[ich][is]<<" "<<rVal<<" "<<rTime<<endl;
803  }
804  cout<<"getDetectorTiming: Event "<<eventseq<<" Gain "<<ig<<" Min "<<minsmpl<<" / "<<rMinVal<<" Max "<<maxsmpl<<" / "<<rMaxVal<<endl;
805  rMinVal -= (ig==HIGH? PEDESTAL*hgDetChannels : PEDESTAL*lgDetChannels);;
806  rMaxVal -= (ig==HIGH? PEDESTAL*hgDetChannels : PEDESTAL*lgDetChannels);;
807  if(abs(rMinVal) > abs(rMaxVal)) {
808  rVal = rMinVal; rTime = minsmpl;
809  } else {
810  rVal = rMaxVal; rTime = maxsmpl;
811  }
812  cout<<"getDetectorTimingT: Signal at "<<rTime<<" Amplitude "<<rVal<<endl;
813 
814  // if(rVal<0.) rVal = 0.;
815  Double_t par0[6];
816  par0[0] = rVal; //3.;
817  par0[1] = max(0.,(Double_t)(rTime-RISETIME));
818  par0[2] = 4.;
819  par0[3] = 1.6;
820  par0[4] = (ig==HIGH? PEDESTAL*hgDetChannels : PEDESTAL*lgDetChannels); // we do not do pedestal subtrastion at this time
821  par0[5] = 0; // slope of the pedestals is initially set to "0"
822  evtShape[ig]->SetParameters(par0);
823  for(int parn=0; parn<6; parn++) evtShape[ig]->SetParLimits(parn, par0Min[parn], par0Max[parn]);
824  evtShape[ig]->SetParLimits(4, par0[4]-PEDESTAL, par0[4]+PEDESTAL);
825  // fitFunc->cd(chan+1); // shapes[chan]->Draw();
826  // gPad->Modified();
827  evtGraph[ig]->Fit(evtShape[ig], "RQWM", "NQ", 0., (Double_t)NSAMPLES);
828  evtShape[ig]->GetParameters(evtfitpar[ig]);
829  cout<<"<getDetectorTiming>: Gain "<<ig<<" Parameters (0) "<<evtfitpar[ig][0]<<" (1) "<<evtfitpar[ig][1]<<" (2) "<<evtfitpar[ig][2]<<" (3) "<<evtfitpar[ig][3]<<" (4) "<<evtfitpar[ig][4]<<" (5) "<<evtfitpar[ig][5]<<endl;
830  evtpedestal[ig] = (evtfitpar[ig][4]+evtfitpar[ig][5]*rTime); //
831  evtSignal[ig]->SetParameters(evtfitpar[ig]);
832  evtSignal[ig]->SetParameter(4, 0.);
833  evtSignal[ig]->SetParameter(5, 0.);
834  Double_t sgp = evtSignal[ig]->GetMaximum(par0Min[1], par0Max[1]);
835  Double_t tmp = evtSignal[ig]->GetMaximumX(par0Min[1], par0Max[1]);
836  Double_t sgn = evtSignal[ig]->GetMinimum(par0Min[1], par0Max[1]);
837  Double_t tmn = evtSignal[ig]->GetMinimumX(par0Min[1], par0Max[1]);
838  evtpeak[ig] = rVal;
839  evttime[ig] = rTime;
840  cout<<"<getDetectorTimingT>: Gain "<<ig<<" Peak0 "<< evtpeak[ig]<<" Time0 "<< evttime[ig]<<endl;
841  if(tmp!=0&&abs(sgp)>abs(sgn)) {
842  evtpeak[ig] = sgp;
843  evttime[ig] = tmp;
844  cout<<"<getDetectorTimingT>: Positive peak selected"<<endl;
845  } else {
846  evtpeak[ig] = sgn;
847  evttime[ig] = tmn;
848  cout<<"<getDetectorTimingT>: Negative peak selected"<<endl;
849  }
850  tChi2[ig] = evtShape[ig]->GetChisquare()/evtShape[ig]->GetNDF();
851  cout<<"<getDetectorTimingT>: Gain "<<ig<<" Peak "<< evtpeak[ig]<<" Time "<< evttime[ig]<<" Chi2 "<<tChi2[ig]<<endl;
852  }
853 
854 }
855 // **************************************************************************
856 
857 void hLabHelper::displaySumEvent(){
858 
859  Int_t ng = (lgDetChannels!=0&&hgDetChannels!=0)? 2 : 1;
860  TCanvas * eventSum = (TCanvas *)(gROOT->FindObject("eventSum"));
861  if(!eventSum) eventSum = new TCanvas("eventSum","sPhenix eventSum display", 300, 300, 800, 500);
862  else eventSum->Clear();
863  eventSum -> Divide(1,ng);
864  for(int ig = 0; ig < ng; ig++){
865  eventSum->cd(ig+1);
866  // for (int is = 0; is < NSAMPLES; is++) evtGraph[ig]->SetPoint(is,(Double_t)is,evtadcsum[ig][is]);
867  // cout<<"<tileDisplay> call to fitTileSignal"<<endl;
868  // eventSum->Divide(1,2);
869  // eventSum->cd(1);
870  // eventSum->SetEditable(kTRUE);
871  evtGraph[ig]->SetMarkerStyle(20);
872  evtGraph[ig]->SetMarkerSize(0.4);
873  eventSum->cd(ig+1);
874  evtGraph[ig]->Draw("acp");
875  // eventSum->cd(2);
876  evtShape[ig]->Draw("same");
877  // evtShape[ig]->Draw();
878  // evtGraph[ig]->Draw("same");
879  eventSum->Update();
880  }
881  eventSum->Update();
882  // eventSum->Clear();
883 }
884 // **************************************************************************
885 // FIXES to data
886 void hLabHelper::dofixes(){
887  // if(runnumber>=1125 && runnumber<1152){
888  // // invert pulse in channel 134 (number 9 in todays parlance
889  // for (int is = 0; is < NSAMPLES; is++) {adc[8][is] -= 2*PEDESTAL; adc[8][is] = abs(adc[8][is]);}
890  // }
891  ;
892 }
893 
894 // **************************************************************************
895 
896 Bool_t hLabHelper::reject(){
897  return false;
898 }
899 // **************************************************************************