Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MyFavoriteMartin.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MyFavoriteMartin.C
1 #include <iostream>
2 #include <cmath>
3 #include <pmonitor.h>
4 #include "MyFavoriteMartin.h"
5 
6 #include "groot.h"
7 #include "ATrace.h"
8 #include "AZigzag.h"
9 #include "FindBlobs.h"
10 #include "FillBlobHist.h"
11 #include "FillHoughHist.h"
12 
13 #include "TH1.h"
14 #include "TH2.h"
15 #include "TH3.h"
16 #include "TNtuple.h"
17 #include "TFile.h"
18 
19 #include "prdfoStream.h"
20 #include "Quiver.h"
21 
22 #include "params.h"
23 
24 //#include <stdlib.h>
25 #include <cstdlib>
26 #include <ctime>
27 #include <map>
28 
29 prdfoStream *os = 0;
30 int init_done = 0;
31 
32 //KD
33 int first_event_nr = 0;
34 //KD
35 
36 using namespace std;
37 
38 // These are the waveforms...
39 // They are the most basic monitor of the SRS system.
40 TH1F *h1a;
41 TH1F *h1b;
42 TH1F *h1c;
43 TH1F *h1d;
44 
45 TH1F *h2a;
46 TH1F *h2b;
47 TH1F *h2c;
48 TH1F *h2d;
49 
50 TH1F *h3a;
51 TH1F *h3b;
52 TH1F *h3c;
53 TH1F *h3d;
54 
55 TH1F *h4a;
56 TH1F *h4b;
57 TH1F *h4c;
58 TH1F *h4d;
59 
60 // These "scope" histograms are the most basic monitor of the DRS4 system
61 TH1F *scope0, *scope1, *scope2, *scope3;
62 
63 // These hnl_D histograms are used to phase the PLL...
64 TH1F *h1a_D;
65 TH1F *h1b_D;
66 TH1F *h1c_D;
67 TH1F *h1d_D;
68 TH1F *h2a_D;
69 TH1F *h2b_D;
70 TH1F *h2c_D;
71 TH1F *h2d_D;
72 TH1F *h3a_D;
73 TH1F *h3b_D;
74 TH1F *h3c_D;
75 TH1F *h3d_D;
76 TH1F *h4a_D;
77 TH1F *h4b_D;
78 TH1F *h4c_D;
79 TH1F *h4d_D;
80 
81 int pinit()
82 {
83  // Here we will book all the histograms, but NONE of the
84  // Ntuples. The reason is that we only want to make the NTuples
85  // when we've made an output file.
86 
87  // Generate the First groot and initialize the Zigzags.
88  groot *Tree = groot::instance();
89 
90  scope0 = new TH1F ( "scope0","waveform Ch 0", 1024, -0.5, 1023.5);
91  scope1 = new TH1F ( "scope1","waveform Ch 1", 1024, -0.5, 1023.5);
92  scope2 = new TH1F ( "scope2","waveform Ch 2", 1024, -0.5, 1023.5);
93  scope3 = new TH1F ( "scope3","waveform Ch 3", 1024, -0.5, 1023.5);
94 
95  scope0->SetLineColor(kYellow+2);
96  scope1->SetLineColor(kBlue);
97  scope2->SetLineColor(kViolet);
98  scope3->SetLineColor(kGreen);
99 
100  h1a = new TH1F( "h1a", "Waveform", 4500, -0.5, 4499.5);
101  h1b = new TH1F( "h1b", "Waveform", 4500, -0.5, 4499.5);
102  h1c = new TH1F( "h1c", "Waveform", 4500, -0.5, 4499.5);
103  h1d = new TH1F( "h1d", "Waveform", 4500, -0.5, 4499.5);
104 
105  h2a = new TH1F( "h2a", "Waveform", 4500, -0.5, 4499.5);
106  h2b = new TH1F( "h2b", "Waveform", 4500, -0.5, 4499.5);
107  h2c = new TH1F( "h2c", "Waveform", 4500, -0.5, 4499.5);
108  h2d = new TH1F( "h2d", "Waveform", 4500, -0.5, 4499.5);
109 
110  h3a = new TH1F( "h3a", "Waveform", 4500, -0.5, 4499.5);
111  h3b = new TH1F( "h3b", "Waveform", 4500, -0.5, 4499.5);
112  h3c = new TH1F( "h3c", "Waveform", 4500, -0.5, 4499.5);
113  h3d = new TH1F( "h3d", "Waveform", 4500, -0.5, 4499.5);
114 
115  h4a = new TH1F( "h4a", "Waveform", 4500, -0.5, 4499.5);
116  h4b = new TH1F( "h4b", "Waveform", 4500, -0.5, 4499.5);
117  h4c = new TH1F( "h4c", "Waveform", 4500, -0.5, 4499.5);
118  h4d = new TH1F( "h4d", "Waveform", 4500, -0.5, 4499.5);
119 
120  h1a_D = new TH1F( "h1a_D", "Waveform", 4001, -0.5, 4000.5);
121  h1b_D = new TH1F( "h1b_D", "Waveform", 4001, -0.5, 4000.5);
122  h1c_D = new TH1F( "h1c_D", "Waveform", 4001, -0.5, 4000.5);
123  h1d_D = new TH1F( "h1d_D", "Waveform", 4001, -0.5, 4000.5);
124  h2a_D = new TH1F( "h2a_D", "Waveform", 4001, -0.5, 4000.5);
125  h2b_D = new TH1F( "h2b_D", "Waveform", 4001, -0.5, 4000.5);
126  h2c_D = new TH1F( "h2c_D", "Waveform", 4001, -0.5, 4000.5);
127  h2d_D = new TH1F( "h2d_D", "Waveform", 4001, -0.5, 4000.5);
128  h3a_D = new TH1F( "h3a_D", "Waveform", 4001, -0.5, 4000.5);
129  h3b_D = new TH1F( "h3b_D", "Waveform", 4001, -0.5, 4000.5);
130  h3c_D = new TH1F( "h3c_D", "Waveform", 4001, -0.5, 4000.5);
131  h3d_D = new TH1F( "h3d_D", "Waveform", 4001, -0.5, 4000.5);
132  h4a_D = new TH1F( "h4a_D", "Waveform", 4001, -0.5, 4000.5);
133  h4b_D = new TH1F( "h4b_D", "Waveform", 4001, -0.5, 4000.5);
134  h4c_D = new TH1F( "h4c_D", "Waveform", 4001, -0.5, 4000.5);
135  h4d_D = new TH1F( "h4d_D", "Waveform", 4001, -0.5, 4000.5);
136 
137  return 0;
138 }
139 
140 
141 int process_event (Event * e)
142 {
143  // Set the run number if event is from new run...
144  int thisRun = e->getRunNumber();
145  if( Quiver::RunNumber != thisRun ) Quiver::RunNumber = thisRun;
146 
147  // Line below is used to put NON-DATA events into output stream
148  if ( os && e->getEvtType() >7 ) os->addEvent(e);
149 
150  // Don't run physics analysis code on NON-DATA events!
151  if (e->getEvtType() > 7) return 0;
152 
153  //Use nr-th event as starting point of process_event
154  // and skip every event before that
155  static int skipentries=0;
156  skipentries++;
157  static bool SetSkipEntry = true;
158  if(SetSkipEntry)
159  {
160  cout << " Skipping " << first_event_nr << " events."
161  << endl;
162  SetSkipEntry = false;
163  }
164  if(skipentries%100==0 && skipentries<first_event_nr)
165  cout << " Skipped " << skipentries << " events." << endl;
166  if(e->getEvtSequence() < first_event_nr && e->getEvtType()==1)
167  return 0;
168 
169  // Here is a data integrity check for the SRS system
170  // Use this for offline analysis at a mature level only.
171  static bool PriorPassed = true;
172  Packet *p = e->getPacket(1010);
173  if (p)
174  {
175  int nhb = p->iValue(0,"NHYBRIDS");
176 
177  // DIagnostic...
178  // cout << "Number of Hybrids: " << nhb << endl;
179  // cout << "Number of Samples: " << endl;
180  // for (int i=0; i<nhb; i++)
181  // {
182  // int nsmp = p->iValue(i,"NSAMPLES");
183  // cout << " " << nsmp;
184  // }
185  // cout << endl;
186 
187  // Check that the data has the right number of hybrids.
188  if (nhb != Nhybrid)
189  {
190  if (PriorPassed)
191  {
192  cout << "Wrong number of hybrids: "
193  << nhb << " found but expected "
194  << Nhybrid << ", abort event."
195  << endl;
196  }
197  else
198  {
199  cout << ".";
200  }
201  delete p;
202  PriorPassed = false;
203  return 0;
204  }
205 
206  // Check that all hybrids have the same data length as the first one.
207  for (int i=1; i<nhb-1; i++)
208  {
209  if (p->iValue(0,"NSAMPLES") != p->iValue(i,"NSAMPLES"))
210  {
211  if (PriorPassed)
212  {
213  cout << "Data length mismatch: "
214  << p->iValue(i,"NSAMPLES")
215  << " out of " << p->iValue(0,"NSAMPLES")
216  << " on hybrid " << i
217  << " abort event."
218  << endl;
219  }
220  else
221  {
222  cout << ".";
223  }
224  delete p;
225  PriorPassed = false;
226  return 0;
227  }
228  }
229  if (!PriorPassed) cout << endl;
230  PriorPassed = true;
231  }
232 
233  // Fancy event counter written to the screen...
234  static int entries=0;
235  entries++;
236  if (entries%10==0)
237  {
238  cout << "Processed " << entries;
239  time_t now = time(0);
240  char* dt = ctime(&now);
241  cout << " at " << dt << endl;
242  }
243 
244  groot *Tree = groot::instance();
245  Tree->ClearTheDetector();
246  Tree->event = e;
247 
248  scope0->Reset();
249  scope1->Reset();
250  scope2->Reset();
251  scope3->Reset();
252 
253  h1a->Reset();
254  h1b->Reset();
255  h1c->Reset();
256  h1d->Reset();
257 
258  h2a->Reset();
259  h2b->Reset();
260  h2c->Reset();
261  h2d->Reset();
262 
263  h3a->Reset();
264  h3b->Reset();
265  h3c->Reset();
266  h3d->Reset();
267 
268  h4a->Reset();
269  h4b->Reset();
270  h4c->Reset();
271  h4d->Reset();
272 
273  // Reading the PSI 4-channel "oscilloscope"...
274  Packet *p1 = e->getPacket(1020);
275  if (p1 && Quiver::AnalyzeScope)
276  {
277  int samples = p1->iValue(0, "SAMPLES");
278  for ( int i = 0; i < samples; i++)
279  {
280 
281  scope0->Fill ( i, -p1->rValue(i,0) );
282  scope1->Fill ( i, -p1->rValue(i,1) );
283  scope2->Fill ( i, -p1->rValue(i,2) );
284  scope3->Fill ( i, -p1->rValue(i,3) );
285 
286  Tree->theTraces[0]->voltage.push_back(-p1->rValue(i,0));
287  Tree->theTraces[1]->voltage.push_back(-p1->rValue(i,1));
288  Tree->theTraces[2]->voltage.push_back(-p1->rValue(i,2));
289  Tree->theTraces[3]->voltage.push_back(-p1->rValue(i,3));
290 
291  }
292  }
293  if (p1) delete p1;
294 
295  // Now monitor and fill SRS to groot.
296  if (p && Quiver::AnalyzeSRS)
297  {
298 
299  int nhb = p->iValue(0,"NHYBRIDS");
300  for (int i=0; i<nhb; i++)
301  {
302  if (p->iValue(i,"NSAMPLES")<24)
303  {
304  cout << "ERROR in Sample Length, "<< p->iValue(i,"NSAMPLES") <<" abort event" << endl;
305  return 0;
306  }
307  }
308 
309  static bool FirstEvent = true;
310  if (FirstEvent)
311  {
312  cout << "Reports for " << nhb << " Hybrids:" << endl;
313  for (int i=0; i<nhb; i++)
314  {
315  cout << i << ": " << p->iValue(i,"NSAMPLES") << endl;
316  }
317  FirstEvent = false;
318  }
319 
320 
321  // Check the hybrid count for sensibility...
322  if (nhb>=1)
323  {
324  int idigi;
325  for (idigi = 0; idigi< 4500; idigi++)
326  {
327  h1a->Fill( idigi, p->iValue(idigi, 0, "RAWSAMPLES"));
328  h1b->Fill( idigi, p->iValue(idigi, 1, "RAWSAMPLES"));
329  h1c->Fill( idigi, p->iValue(idigi, 2, "RAWSAMPLES"));
330  h1d->Fill( idigi, p->iValue(idigi, 3, "RAWSAMPLES"));
331 
332  h2a->Fill( idigi, p->iValue(idigi, 4, "RAWSAMPLES"));
333  h2b->Fill( idigi, p->iValue(idigi, 5, "RAWSAMPLES"));
334  h2c->Fill( idigi, p->iValue(idigi, 6, "RAWSAMPLES"));
335  h2d->Fill( idigi, p->iValue(idigi, 7, "RAWSAMPLES"));
336 
337  h3a->Fill( idigi, p->iValue(idigi, 8, "RAWSAMPLES"));
338  h3b->Fill( idigi, p->iValue(idigi, 9, "RAWSAMPLES"));
339  h3c->Fill( idigi, p->iValue(idigi,10, "RAWSAMPLES"));
340  h3d->Fill( idigi, p->iValue(idigi,11, "RAWSAMPLES"));
341 
342  h4a->Fill( idigi, p->iValue(idigi, 12, "RAWSAMPLES"));
343  h4b->Fill( idigi, p->iValue(idigi, 13, "RAWSAMPLES"));
344  h4c->Fill( idigi, p->iValue(idigi, 14, "RAWSAMPLES"));
345  h4d->Fill( idigi, p->iValue(idigi, 15, "RAWSAMPLES"));
346 
347  if (idigi>500 && idigi<3500)
348  {
349  h1a_D->Fill( p->iValue(idigi, 0, "RAWSAMPLES"));
350  h1b_D->Fill( p->iValue(idigi, 1, "RAWSAMPLES"));
351  h1c_D->Fill( p->iValue(idigi, 2, "RAWSAMPLES"));
352  h1d_D->Fill( p->iValue(idigi, 3, "RAWSAMPLES"));
353  h2a_D->Fill( p->iValue(idigi, 4, "RAWSAMPLES"));
354  h2b_D->Fill( p->iValue(idigi, 5, "RAWSAMPLES"));
355  h2c_D->Fill( p->iValue(idigi, 6, "RAWSAMPLES"));
356  h2d_D->Fill( p->iValue(idigi, 7, "RAWSAMPLES"));
357  h3a_D->Fill( p->iValue(idigi, 8, "RAWSAMPLES"));
358  h3b_D->Fill( p->iValue(idigi, 9, "RAWSAMPLES"));
359  h3c_D->Fill( p->iValue(idigi,10, "RAWSAMPLES"));
360  h3d_D->Fill( p->iValue(idigi,11, "RAWSAMPLES"));
361  h4a_D->Fill( p->iValue(idigi,12, "RAWSAMPLES"));
362  h4b_D->Fill( p->iValue(idigi,13, "RAWSAMPLES"));
363  h4c_D->Fill( p->iValue(idigi,14, "RAWSAMPLES"));
364  h4d_D->Fill( p->iValue(idigi,15, "RAWSAMPLES"));
365  }
366  }
367 
368  // Read the tracker raw data into the raw array...
369  int rindex=0;
370  for (int JINX =0; JINX<nhb; JINX++)
371  {
372  for ( int j = 0; j< CHhybrid; j++)
373  {
374  for (int i=0; i<min(p->iValue(JINX,"NSAMPLES"),24); i++) // limit to 24 samples...
375  {
376  // For now we'll use a trivial routine to analyze waveform.
377  // We'll just sum across the samples...
378  AZigzag::Raw[rindex].push_back(p->iValue(j,i,JINX));
379  }
380  rindex++;
381  }
382  }
383  AZigzag::DetermineCommonMode();
384  AZigzag::ApplyCalibration();
385  for (int i=0; i<Tree->theZigzags.size(); i++)
386  {
387  Tree->theZigzags[i]->DetermineQ();
388  }
389  FindBlobs();
390  FillBlobHist();
391  //FillHoughHist();
392  }
393  else
394  {
395  cout << "ERROR " << nhb << " Hybrids Found:" << endl;
396  for (int i=0; i<nhb; i++)
397  {
398  cout << i << ": " << p->iValue(i,"NSAMPLES") << endl;
399  }
400  }
401  }
402 
403 
404  // Here select whether to write the event to the output stream.
405  static int nSaved=0;
406  if ( os && Tree->theTracks.size() >= 1 )
407  {
408  nSaved++;
409  os->addEvent(e);
410  if (nSaved%100==0)
411  {
412  cout << "Saved " << nSaved;
413  cout << " percentage "
414  << 100.0*(double)nSaved/(double)entries;
415  cout << endl;
416  }
417  }
418  if (p) delete p;
419 
420  return 0;
421 }
422 
423 
424 
425 int prdfopen ( const char *filename)
426 {
427  if (os)
428  {
429  cout << "file is already open " << endl;
430  return -1;
431  }
432 
433  os = new prdfoStream (filename);
434  if ( os->is_defunct() )
435  {
436  delete os;
437  os = 0;
438  return -1;
439  }
440 
441  return 0;
442 }
443 
444 int prdfclose ()
445 {
446  if ( !os)
447  {
448  cout << "no file open " << endl;
449  return -1;
450  }
451  delete os;
452  os = 0;
453  return 0;
454 }
455 
456 
457 
458 //Function that allows to start prun() from the nr-th event
459 int getFirstEventNr() { return first_event_nr; }
460 
461 void setFirstEventNr(const int nr) { first_event_nr = nr; }
462 
463 // Binary number to GrayCode (courtesy of Martin)
464 int b2G(int num) { return (num >> 1) ^ num; }