[u/mrichter/AliRoot.git] / EMCAL / SMcalib / QA / LEDRef_evtdis.C

///////////////////////////////////////////////
// Simple event display for channel spectra
// Input file should be raw DATE root tree file
///////////////////////////////////////////////

// for the readout: 1 FEE reads out LED ref. info for the whole SuperModule
// 1 CSP per StripModule
const int NSTRIPS = 24; // number of StripModules in SuperModule
const int TOTCHAN = NSTRIPS * 2;	// *2 since we have high gain and low gain
// we group the calibrations in 3 sets of 8 strips
const int NSETS = 3; 
const int NSTRIPS_IN_SET = 8;

// gamma2 fit function
double fitfun(double *x, double *par) {
  double Amp	= par[0];
  double Tmax	= par[1];
  double Tau	= par[2];
  double Ped	= par[3];
  double gammaN   = par[4];
  double t = 0;
  if(Tau) t = (x[0] - Tmax + Tau)/Tau;
  if(t<0) t = 0;
  
  // Since we have now set gammaN to 2, we replace the pow(t, 2) call in
  // double f = Amp * pow(t,gammaN) * exp(gammaN*(1-t)) + Ped;
  // with just t*t
  double f = Amp * t*t * exp(gammaN*(1-t)) + Ped;
  return f;
}

// main method
void LEDRef_evtdis(const int runno = 615,
		   const int gainv = 0,  /*0=low, 1=high*/
		   const int evtnum= -10,
		   int ymax=1023, // set the scale of plots
		   const int delay = 1)  // -1=no delay, wait for input, X>=0 => sleep aprox. X sec. after making plot
{
  // set ranges to plot
  const int strip_f = 0; // first
  const int strip_l = NSTRIPS - 1;
  
  const int nsamples = 65; // number of ADC time samples per channel and event
  
  const int saveplot = 0;
  const int numbering      = 1; // 0: no numbering, 1: nubering on each plot
  const int dofit = 0; // 0: no fit, 1: try to fit the spectra 
  const int debug    = 0;
  const float gammaN = 2;
  // end of setup    
  
  // Assume we are just interested in the 1st segment, _0.root below for fname*
  Char_t fname[256];
  sprintf(fname, "/local/data/Run_%09d.Seq_1A.Stream_0.root",runno);
  cout << "TOTCHAN " << TOTCHAN << endl;

  // set up a raw reader of the data
  AliRawReader *rawReader = NULL;
  rawReader = new AliRawReaderRoot(fname);
  AliCaloRawStream *in = NULL; 
  in = new AliCaloRawStream(rawReader,"EMCAL");

  // set up histograms
  TH1F *hfit[TOTCHAN];
  TF1 *f1[TOTCHAN];
  char ch_label[TOTCHAN][100];
  char buff1[100];
  char name[80];
  for(int i=0; i<TOTCHAN; i++) {
    sprintf(buff1,"hfit_%d",i);
    hfit[i] = new TH1F(buff1,"hfit", nsamples , -0.5, nsamples - 0.5);
    hfit[i]->SetDirectory(0);
    sprintf(name,"f1_%d",i);
    f1[i] = new TF1(name,fitfun,0,70,5);
    f1[i]->SetLineWidth(2);
    f1[i]->SetLineColor(2);

    //	int idx = istrip + NSTRIPS * gain; // encoding used later
    int gain = i / (NSTRIPS);
    int istrip = i % NSTRIPS;
    sprintf(ch_label[i], "Strip%02d", istrip);
  }
  
  TCanvas *cc1 = new TCanvas("cc1","3 columns of 8 strips each",600,800);
  int numcol = NSETS;
  int numrow = NSTRIPS_IN_SET;
  cc1->Divide(numcol, numrow);
  
  TText *t = new TText;
  t->SetTextSize(0.17);
  int clr[2] = {4,2}; // colors
  
  // figure out which events we should look at
  int firstevent = evtnum;
  int lastevent = evtnum;
  if (evtnum < 0) { // get a bunch of events
    firstevent = 0;
    lastevent = - evtnum;
  }
  if (evtnum == 0) { // get all events
    firstevent = 0;
    lastevent = 1000000;
  }
  
  Int_t iev =0;
  AliRawEventHeaderBase *aliHeader=NULL;    
  while ( rawReader->NextEvent() && iev < firstevent) {
    aliHeader = (AliRawEventHeaderBase*) rawReader->GetEventHeader();
    iev++;
  }
  
  // loop over selected events
  while ( rawReader->NextEvent() && iev <= lastevent) {
    aliHeader = (AliRawEventHeaderBase*) rawReader->GetEventHeader();
    int runNumber = aliHeader->Get("RunNb"); 
    
    cout << "Found run number " << runNumber << endl;
    
    // reset histograms
    for(int i=0; i<TOTCHAN; i++) {
      hfit[i]->Reset();
    }
    
    // get events (the "1" ensures that we actually select all events for now)
    if ( 1 || aliHeader->Get("Type") == AliRawEventHeaderBase::kPhysicsEvent ) {
      const UInt_t * evtId = aliHeader->GetP("Id");
      int evno_raw = (int) evtId[0];
      int timestamp = aliHeader->Get("Timestamp");
      
      cout << " evno " << evno_raw
	   << " size " << aliHeader->GetEventSize()
	   << " type " << aliHeader->Get("Type")
	   << " type name " << aliHeader->GetTypeName()
	   << " timestamp " << timestamp
	   << endl;
      
      /// process_event stream
      while ( in->Next() ) {
	
	int strip = in->GetColumn();
	int gain = in->GetRow();
	
	if (in->IsLEDMonData()) {
	  
	  int idx = strip + NSTRIPS*gain;
	  //cout << "hist idx " << idx << endl;
	  
	  if (idx < 0 || idx > TOTCHAN) { 
	    cout << "Hist idx out of range: " << idx << endl;
	  }
	  else { // reasonable range of idx
	    hfit[idx]->SetBinContent(in->GetTime(), in->GetSignal());
	  }

	} // LED Ref data only

      } // Raw data read
    
      // Next: let's actually plot the data..
      for (Int_t strip = strip_f; strip <= strip_l; strip++) {
	
	int idx = strip + NSTRIPS*gainv;
	
	// which set/column does the strip belong in
	int iset = strip / NSTRIPS_IN_SET; 	  
	int within_set = strip % NSTRIPS_IN_SET; 	  
	// on which pad should we plot it?
	int pad_id = (NSTRIPS_IN_SET-1-within_set)*NSETS + iset + 1;
	
	cout << "strip " << strip 
	     << ". set="<< iset << ", within_set=" << within_set
	     << ", pad=" << pad_id << endl;
	cc1->cd(pad_id);
	hfit[idx]->SetTitle("");
	hfit[idx]->SetFillColor(5);
	hfit[idx]->SetMaximum(ymax);
	hfit[idx]->SetMinimum(0);
	// we may or may not decide to fit the data
	if (dofit) {
	  f1[i]->SetParameter(0, 0); // initial guess; zero amplitude :=)
	  hfit[idx]->Fit(f1[i]);
	}
	hfit[idx]->Draw();
	if( numbering ) {
	  t->SetTextColor(clr[gainv]);
	  t->DrawTextNDC(0.65,0.65,ch_label[idx]);
	}
      }

      // add some extra text on the canvas
      // print a box showing run #, evt #, and timestamp
      cc1->cd();
      // first draw transparent pad
      TPad *trans = new TPad("trans","",0,0,1,1);
      trans->SetFillStyle(4000);
      trans->Draw();
      trans->cd();
      // then draw text
      TPaveText *label = new TPaveText(.2,.11,.8,.14,"NDC"); 
      //  label->Clear();
      label->SetBorderSize(1);
      label->SetFillColor(0);
      label->SetLineColor(clr[gainv]);
      label->SetTextColor(clr[gainv]);
      //label->SetFillStyle(0);
      TDatime d;
      d.Set(timestamp);
      sprintf(name,"Run %d, Event %d, Hist Max %d, %s",runno,iev,ymax,d.AsString());
      label->AddText(name);
      label->Draw();
      cc1->Update();
      cout << "Done" << endl;
      
      // some shenanigans to hold the plotting, if requested
      if (firstevent != lastevent) {
	if (delay == -1) {
	  // wait for character input before proceeding
	  cout << " enter y to proceed " << endl;
	  char dummy[2];
	  cin >> dummy;
	  cout << " read " << dummy << endl;
	  if (strcmp(dummy, "y")==0) {
	    cout << " ok, continuing with event " << iev+1 << endl;
	  }
	  else {
	    cout << " ok, exiting " << endl;
	    //exit(1);
	  }
	}
	else {
	  cout << "Sleeping for " << delay * 500 << endl;
	  gSystem->Sleep(delay * 500);
	}
      }

      // save plot, if setup/requested to do so
      char plotname[100];
      if (saveplot==1) {
	sprintf(plotname,"Run_%d_LEDRef_Ev%d_Gain%d_MaxHist%d.gif",
		runno,iev,gainv,ymax);  
	cout <<"SAVING plot:"<< plotname << endl;
	cc1->SaveAs(plotname);
      }

    } // event selection

  iev ++;
  } // event loop

}
Commit	Line	Data
d63882b1	1	///////////////////////////////////////////////
	2	// Simple event display for channel spectra
	3	// Input file should be raw DATE root tree file
	4	///////////////////////////////////////////////
	5
	6	// for the readout: 1 FEE reads out LED ref. info for the whole SuperModule
	7	// 1 CSP per StripModule
	8	const int NSTRIPS = 24; // number of StripModules in SuperModule
	9	const int TOTCHAN = NSTRIPS * 2; // *2 since we have high gain and low gain
	10	// we group the calibrations in 3 sets of 8 strips
	11	const int NSETS = 3;
	12	const int NSTRIPS_IN_SET = 8;
	13
	14	// gamma2 fit function
	15	double fitfun(double x, double par) {
	16	double Amp = par[0];
	17	double Tmax = par[1];
	18	double Tau = par[2];
	19	double Ped = par[3];
	20	double gammaN = par[4];
	21	double t = 0;
	22	if(Tau) t = (x[0] - Tmax + Tau)/Tau;
	23	if(t<0) t = 0;
	24
	25	// Since we have now set gammaN to 2, we replace the pow(t, 2) call in
	26	// double f = Amp * pow(t,gammaN) * exp(gammaN*(1-t)) + Ped;
	27	// with just t*t
	28	double f = Amp * tt exp(gammaN*(1-t)) + Ped;
	29	return f;
	30	}
	31
	32	// main method
	33	void LEDRef_evtdis(const int runno = 615,
	34	const int gainv = 0, /0=low, 1=high/
	35	const int evtnum= -10,
	36	int ymax=1023, // set the scale of plots
	37	const int delay = 1) // -1=no delay, wait for input, X>=0 => sleep aprox. X sec. after making plot
	38	{
	39	// set ranges to plot
	40	const int strip_f = 0; // first
	41	const int strip_l = NSTRIPS - 1;
	42
	43	const int nsamples = 65; // number of ADC time samples per channel and event
	44
	45	const int saveplot = 0;
	46	const int numbering = 1; // 0: no numbering, 1: nubering on each plot
	47	const int dofit = 0; // 0: no fit, 1: try to fit the spectra
	48	const int debug = 0;
	49	const float gammaN = 2;
	50	// end of setup
	51
	52	// Assume we are just interested in the 1st segment, _0.root below for fname*
	53	Char_t fname[256];
	54	sprintf(fname, "/local/data/Run_%09d.Seq_1A.Stream_0.root",runno);
	55	cout << "TOTCHAN " << TOTCHAN << endl;
	56
	57	// set up a raw reader of the data
	58	AliRawReader *rawReader = NULL;
	59	rawReader = new AliRawReaderRoot(fname);
	60	AliCaloRawStream *in = NULL;
	61	in = new AliCaloRawStream(rawReader,"EMCAL");
	62
	63	// set up histograms
	64	TH1F *hfit[TOTCHAN];
65	TF1 *f1[TOTCHAN];
66	char ch_label[TOTCHAN][100];
67	char buff1[100];
68	char name[80];
69	for(int i=0; i<TOTCHAN; i++) {
70	sprintf(buff1,"hfit_%d",i);
71	hfit[i] = new TH1F(buff1,"hfit", nsamples , -0.5, nsamples - 0.5);
72	hfit[i]->SetDirectory(0);
73	sprintf(name,"f1_%d",i);
74	f1[i] = new TF1(name,fitfun,0,70,5);
75	f1[i]->SetLineWidth(2);
76	f1[i]->SetLineColor(2);
77
78	// int idx = istrip + NSTRIPS * gain; // encoding used later
79	int gain = i / (NSTRIPS);
80	int istrip = i % NSTRIPS;
81	sprintf(ch_label[i], "Strip%02d", istrip);
82	}
83
84	TCanvas *cc1 = new TCanvas("cc1","3 columns of 8 strips each",600,800);
85	int numcol = NSETS;
86	int numrow = NSTRIPS_IN_SET;
87	cc1->Divide(numcol, numrow);
88
89	TText *t = new TText;
90	t->SetTextSize(0.17);
91	int clr[2] = {4,2}; // colors
92
93	// figure out which events we should look at
94	int firstevent = evtnum;
95	int lastevent = evtnum;
96	if (evtnum < 0) { // get a bunch of events
97	firstevent = 0;
98	lastevent = - evtnum;
99	}
100	if (evtnum == 0) { // get all events
101	firstevent = 0;
102	lastevent = 1000000;
103	}
104
105	Int_t iev =0;
106	AliRawEventHeaderBase *aliHeader=NULL;
107	while ( rawReader->NextEvent() && iev < firstevent) {
108	aliHeader = (AliRawEventHeaderBase*) rawReader->GetEventHeader();
109	iev++;
110	}
111
112	// loop over selected events
113	while ( rawReader->NextEvent() && iev <= lastevent) {
114	aliHeader = (AliRawEventHeaderBase*) rawReader->GetEventHeader();
115	int runNumber = aliHeader->Get("RunNb");
116
117	cout << "Found run number " << runNumber << endl;
118
119	// reset histograms
120	for(int i=0; i<TOTCHAN; i++) {
121	hfit[i]->Reset();
122	}
123
124	// get events (the "1" ensures that we actually select all events for now)
125	if ( 1 \|\| aliHeader->Get("Type") == AliRawEventHeaderBase::kPhysicsEvent ) {
126	const UInt_t * evtId = aliHeader->GetP("Id");
127	int evno_raw = (int) evtId[0];
128	int timestamp = aliHeader->Get("Timestamp");
129
130	cout << " evno " << evno_raw
131	<< " size " << aliHeader->GetEventSize()
132	<< " type " << aliHeader->Get("Type")
133	<< " type name " << aliHeader->GetTypeName()
134	<< " timestamp " << timestamp
135	<< endl;
136
137	/// process_event stream
138	while ( in->Next() ) {
139
140	int strip = in->GetColumn();
141	int gain = in->GetRow();
142
143	if (in->IsLEDMonData()) {
144
145	int idx = strip + NSTRIPS*gain;
146	//cout << "hist idx " << idx << endl;
147
148	if (idx < 0 \|\| idx > TOTCHAN) {
149	cout << "Hist idx out of range: " << idx << endl;
150	}
151	else { // reasonable range of idx
152	hfit[idx]->SetBinContent(in->GetTime(), in->GetSignal());
153	}
154
155	} // LED Ref data only
156
157	} // Raw data read
158
159	// Next: let's actually plot the data..
160	for (Int_t strip = strip_f; strip <= strip_l; strip++) {
161
162	int idx = strip + NSTRIPS*gainv;
163
164	// which set/column does the strip belong in
165	int iset = strip / NSTRIPS_IN_SET;
166	int within_set = strip % NSTRIPS_IN_SET;
167	// on which pad should we plot it?
168	int pad_id = (NSTRIPS_IN_SET-1-within_set)*NSETS + iset + 1;
169
170	cout << "strip " << strip
171	<< ". set="<< iset << ", within_set=" << within_set
172	<< ", pad=" << pad_id << endl;
173	cc1->cd(pad_id);
174	hfit[idx]->SetTitle("");
175	hfit[idx]->SetFillColor(5);
176	hfit[idx]->SetMaximum(ymax);
177	hfit[idx]->SetMinimum(0);
178	// we may or may not decide to fit the data
179	if (dofit) {
180	f1[i]->SetParameter(0, 0); // initial guess; zero amplitude :=)
181	hfit[idx]->Fit(f1[i]);
182	}
183	hfit[idx]->Draw();
184	if( numbering ) {
185	t->SetTextColor(clr[gainv]);
186	t->DrawTextNDC(0.65,0.65,ch_label[idx]);
187	}
188	}
189
190	// add some extra text on the canvas
191	// print a box showing run #, evt #, and timestamp
192	cc1->cd();
193	// first draw transparent pad
194	TPad *trans = new TPad("trans","",0,0,1,1);
195	trans->SetFillStyle(4000);
196	trans->Draw();
197	trans->cd();
198	// then draw text
199	TPaveText *label = new TPaveText(.2,.11,.8,.14,"NDC");
200	// label->Clear();
201	label->SetBorderSize(1);
202	label->SetFillColor(0);
203	label->SetLineColor(clr[gainv]);
204	label->SetTextColor(clr[gainv]);
205	//label->SetFillStyle(0);
206	TDatime d;
207	d.Set(timestamp);
208	sprintf(name,"Run %d, Event %d, Hist Max %d, %s",runno,iev,ymax,d.AsString());
209	label->AddText(name);
210	label->Draw();
211	cc1->Update();
212	cout << "Done" << endl;
213
214	// some shenanigans to hold the plotting, if requested
215	if (firstevent != lastevent) {
216	if (delay == -1) {
217	// wait for character input before proceeding
218	cout << " enter y to proceed " << endl;
219	char dummy[2];
220	cin >> dummy;
221	cout << " read " << dummy << endl;
222	if (strcmp(dummy, "y")==0) {
223	cout << " ok, continuing with event " << iev+1 << endl;
224	}
225	else {
226	cout << " ok, exiting " << endl;
227	//exit(1);
228	}
229	}
230	else {
231	cout << "Sleeping for " << delay * 500 << endl;
232	gSystem->Sleep(delay * 500);
233	}
234	}
235
236	// save plot, if setup/requested to do so
237	char plotname[100];
238	if (saveplot==1) {
239	sprintf(plotname,"Run_%d_LEDRef_Ev%d_Gain%d_MaxHist%d.gif",
240	runno,iev,gainv,ymax);
241	cout <<"SAVING plot:"<< plotname << endl;
242	cc1->SaveAs(plotname);
243	}
244
245	} // event selection
246
247	iev ++;
248	} // event loop
249
250	}