[u/mrichter/AliRoot.git] / PHOS / macros / RawProcessing / PedAmp.C

#if !defined(__CINT__) || defined(__MAKECINT__)
#include "AliRawReaderRoot.h"
#include "AliCaloRawStream.h"
#include "TString.h"
#include "TH2F.h"
#include "TH1S.h"
#include "TFile.h"
#include "TStopwatch.h"
#include "iostream.h"

#endif

/* $Id$ */

void PedAmp(const char *rawFile, const Int_t debug=0)
{
  // Read raw data, decode it to samples,
  // calculate pedestals from presamples, 
  // evaluate the signal amplitude as a maximum sample, 
  // and fill histograms with pedestals and amplitudes
  // This script should be compiled to speed up the data processing:
  // .L PedAmp.C++
  //___
  // Yuri Kharlov. 6 September 2007

  TStopwatch stopwatch;
  stopwatch.Start();
  
  AliRawReaderRoot* rf = new AliRawReaderRoot(rawFile);
  AliCaloRawStream in(rf,"PHOS");
  in.SetOldRCUFormat(kTRUE);
  const Int_t nPresample = 10;

  TString baseNamePed ="hPed";
  TString baseTitlePed="Ped in cell (";
  TString baseNameLed ="hLED";
  TString baseTitleLed="LED in cell (";
  const char* sgain[2]={"high","low"};

  const Int_t gainMax=2,modMax=5,colMax=56,rowMax=64;
  TH1F *hPed[2][56][64];
  TH1F *hLed[2][56][64];
  for (Int_t gain=0; gain<gainMax; gain++) {
    for (Int_t mod=0; mod<modMax; mod++) {
      for (Int_t col=0; col<colMax; col++) {
	for (Int_t row=0; row<rowMax; row++) {
	  hPed[gain][col][row] = 0;
	  hLed[gain][col][row] = 0;
	}
      }
    }
  }
  TH1F *hPedHiMean1 = new TH1F("hPedHiMean1","Mean pedestals, high gain" ,100,0.,100.);
  TH1F *hPedHiRMS1  = new TH1F("hPedHiRMS1" ,"RMS pedestals, high gain"  ,100,0.,50.);
  TH1F *hPedLoMean1 = new TH1F("hPedLoMean1","Mean pedestals, low gain"  ,100,0.,100.);
  TH1F *hPedLoRMS1  = new TH1F("hPedLoRMS1" ,"RMS pedestals, low gain"   ,100,0.,50.);
  hPedHiMean1->Sumw2();
  hPedHiRMS1 ->Sumw2();
  hPedLoMean1->Sumw2();
  hPedLoRMS1 ->Sumw2();

  TH2F *hPedHiMean  = new TH2F("hPedHiMean","Mean pedestals, high gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hPedHiRMS   = new TH2F("hPedHiRMS" ,"R.M.S. of pedestals, high gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hPedHiNum   = new TH2F("hPedHiNum" ,"Number of pedestals, high gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hPedLoMean  = new TH2F("hPedLoMean","Mean pedestals, low gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hPedLoRMS   = new TH2F("hPedLoRMS" ,"R.M.S. of pedestals, low gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hPedLoNum   = new TH2F("hPedLoNum" ,"Number of pedestals, low gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  
  TH2F *hLedHiMean  = new TH2F("hLedHiMean","Mean LED, high gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hLedHiRMS   = new TH2F("hLedHiRMS" ,"R.M.S. of LED, high gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hLedHiNum   = new TH2F("hLedHiNum" ,"Number of LED, high gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hLedLoMean  = new TH2F("hLedLoMean","Mean LED, low gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hLedLoRMS   = new TH2F("hLedLoRMS" ,"R.M.S. of LED, low gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);
  TH2F *hLedLoNum   = new TH2F("hLedLoNum" ,"Number of LED, low gain",
			       rowMax,0.,rowMax, colMax,0.,colMax);

  Int_t iEvent=0;
  Int_t iBin=0;
  Int_t gain=-1,mod=-1,col=-1,row=-1;
  Int_t runNum = 0;
  Double_t maxAmp = 0;
//   Double_t meanPed,rmsPed;
  Float_t pedMean = 0;
  Int_t   nPed = 0;
  while (rf->NextEvent()) {
    cout << "\n\n\tEvent "<< iEvent << " of type " << rf->GetType() << endl;
    runNum = rf->GetRunNumber();
    while ( in.Next() ) {
      if (!in.IsNewHWAddress() && debug > 1) {
	cout << " time=" << in.GetTime() << "/"<< in.GetTimeLength()
	     << " amp=" <<  in.GetSignal() <<endl;
      }
      if (in.IsNewHWAddress()) {
	iBin   = 0;
	if (gain!=-1 && col!=-1 && row!=-1) {
	  if (nPed < 1) continue;
	  maxAmp -= (Double_t)(pedMean/nPed); // "pedestal subtraction"
	  if (maxAmp > 5) hLed[gain][col][row]->Fill(maxAmp);
	}
	maxAmp  = 0;
	pedMean = 0;
	nPed    = 0;
	gain = in.IsLowGain(); // 0-high, 1-low
	mod  = in.GetModule();
	col  = in.GetColumn();
	row  = in.GetRow();
	if (debug > 0) {
	  cout << "\n\tNew HW address: "<<in.GetHWAddress()
	       << ", gain="<< gain
	       << ", module="<<mod
	       << ", xz=("<<col<<","<<row<<")"<< endl;
	}

	// Create a new histo for the new cell, if it does not exist yet
	if (!hPed[gain][col][row]) {
	  TString name  = baseNamePed;
	  TString title = baseTitlePed;
	  name +="_g"; name +=gain;
	  name +="_m"; name +=mod;
	  name +="_z"; name +=col;
	  name +="_x"; name +=row;

	  title +=mod; title +=",";
	  title +=col; title +=",";
	  title +=row; title +="), ";
	  title +=sgain[gain]; title +=" gain";

	  hPed[gain][col][row] = new TH1F(name,title,100,0.,100.);
	  hPed[gain][col][row]->Sumw2();
	  hPed[gain][col][row]->SetMarkerStyle(20);
	  hPed[gain][col][row]->SetOption("eph");
	}
	if (!hLed[gain][col][row]) {
	  TString name  = baseNameLed;
	  TString title = baseTitleLed;
	  name +="_g"; name +=gain;
	  name +="_m"; name +=mod;
	  name +="_z"; name +=col;
	  name +="_x"; name +=row;

	  title +=mod; title +=",";
	  title +=col; title +=",";
	  title +=row; title +="), ";
	  title +=sgain[gain]; title +=" gain";

	  hLed[gain][col][row] = new TH1F(name,title,1023,0.,1023.);
	  hLed[gain][col][row]->Sumw2();
	  hLed[gain][col][row]->SetMarkerStyle(20);
	  hLed[gain][col][row]->SetOption("eph");
	}
      }

      iBin++;
      // Assume that first nPresample are pedestals
      if (in.GetTimeLength()-iBin < nPresample) {
	hPed[gain][col][row]->Fill(in.GetSignal());
	pedMean += in.GetSignal();
	nPed++;
      }
      if (in.GetSignal() > maxAmp) maxAmp = in.GetSignal();
    } // end of next signal
    iEvent++;
  } // end of next event

  // Fill 2-dim histograms for mean, rms and n pedestals

  gain = 0;
  for (Int_t col=0; col<colMax; col++) {
    for (Int_t row=0; row<rowMax; row++) {
      if (hPed[gain][col][row] != 0) {
	hPedHiMean1->Fill( hPed[gain][col][row]->GetMean());
	hPedHiRMS1 ->Fill( hPed[gain][col][row]->GetRMS() );
	hPedHiMean ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetMean()    );
	hPedHiRMS  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetRMS()     );
	hPedHiNum  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetEntries() );
      }
      if (hLed[gain][col][row] != 0) {
	hLedHiMean ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetMean()    );
	hLedHiRMS  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetRMS()     );
	hLedHiNum  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetEntries() );
      }
    }
  }
  gain = 1;
  for (Int_t col=0; col<colMax; col++) {
    for (Int_t row=0; row<rowMax; row++) {
      if (hPed[gain][col][row] != 0) {
	hPedLoMean1->Fill( hPed[gain][col][row]->GetMean());
	hPedLoRMS1 ->Fill( hPed[gain][col][row]->GetRMS() );
	hPedLoMean ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetMean()    );
	hPedLoRMS  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetRMS()     );
	hPedLoNum  ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetEntries() );
      }
      if (hLed[gain][col][row] != 0) {
	hLedLoMean ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetMean()    );
	hLedLoRMS  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetRMS()     );
	hLedLoNum  ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetEntries() );
      }
    }
  }

  // Write existing histograms to a root file

  TString fileName = "ped";
  fileName += runNum;
  fileName += ".root";
  TFile *file = new TFile(fileName,"RECREATE");

  for (Int_t gain=0; gain<gainMax; gain++) {
    for (Int_t col=0; col<colMax; col++) {
      for (Int_t row=0; row<rowMax; row++) {
	if (hPed[gain][col][row] != 0)
	  hPed[gain][col][row]->Write();
	if (hLed[gain][col][row] != 0)
	  hLed[gain][col][row]->Write();
      }
    }
  }
  hPedHiMean1->Write();
  hPedHiRMS1 ->Write();
  hPedLoMean1->Write();
  hPedLoRMS1 ->Write();
  hPedHiMean ->Write();
  hPedHiRMS  ->Write();
  hPedHiNum  ->Write();
  hPedLoMean ->Write();
  hPedLoRMS  ->Write();
  hPedLoNum  ->Write();
  
  hLedHiMean ->Write();
  hLedHiRMS  ->Write();
  hLedHiNum  ->Write();
  hLedLoMean ->Write();
  hLedLoRMS  ->Write();
  hLedLoNum  ->Write();
  
  file->Close();
  stopwatch.Print();
}
Commit	Line	Data
03b71a7c	1	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	2	#include "AliRawReaderRoot.h"
	3	#include "AliCaloRawStream.h"
	4	#include "TString.h"
	5	#include "TH2F.h"
	6	#include "TH1S.h"
	7	#include "TFile.h"
	8	#include "TStopwatch.h"
	9	#include "iostream.h"
	10
	11	#endif
	12
	13	/* $Id$ */
	14
	15	void PedAmp(const char *rawFile, const Int_t debug=0)
	16	{
	17	// Read raw data, decode it to samples,
	18	// calculate pedestals from presamples,
	19	// evaluate the signal amplitude as a maximum sample,
	20	// and fill histograms with pedestals and amplitudes
	21	// This script should be compiled to speed up the data processing:
	22	// .L PedAmp.C++
	23	//___
	24	// Yuri Kharlov. 6 September 2007
	25
	26	TStopwatch stopwatch;
	27	stopwatch.Start();
	28
	29	AliRawReaderRoot* rf = new AliRawReaderRoot(rawFile);
	30	AliCaloRawStream in(rf,"PHOS");
	31	in.SetOldRCUFormat(kTRUE);
	32	const Int_t nPresample = 10;
	33
	34	TString baseNamePed ="hPed";
	35	TString baseTitlePed="Ped in cell (";
	36	TString baseNameLed ="hLED";
	37	TString baseTitleLed="LED in cell (";
	38	const char* sgain[2]={"high","low"};
	39
	40	const Int_t gainMax=2,modMax=5,colMax=56,rowMax=64;
	41	TH1F *hPed[2][56][64];
	42	TH1F *hLed[2][56][64];
	43	for (Int_t gain=0; gain<gainMax; gain++) {
	44	for (Int_t mod=0; mod<modMax; mod++) {
	45	for (Int_t col=0; col<colMax; col++) {
	46	for (Int_t row=0; row<rowMax; row++) {
	47	hPed[gain][col][row] = 0;
	48	hLed[gain][col][row] = 0;
	49	}
	50	}
	51	}
	52	}
	53	TH1F *hPedHiMean1 = new TH1F("hPedHiMean1","Mean pedestals, high gain" ,100,0.,100.);
	54	TH1F *hPedHiRMS1 = new TH1F("hPedHiRMS1" ,"RMS pedestals, high gain" ,100,0.,50.);
	55	TH1F *hPedLoMean1 = new TH1F("hPedLoMean1","Mean pedestals, low gain" ,100,0.,100.);
	56	TH1F *hPedLoRMS1 = new TH1F("hPedLoRMS1" ,"RMS pedestals, low gain" ,100,0.,50.);
	57	hPedHiMean1->Sumw2();
	58	hPedHiRMS1 ->Sumw2();
	59	hPedLoMean1->Sumw2();
	60	hPedLoRMS1 ->Sumw2();
	61
	62	TH2F *hPedHiMean = new TH2F("hPedHiMean","Mean pedestals, high gain",
	63	rowMax,0.,rowMax, colMax,0.,colMax);
	64	TH2F *hPedHiRMS = new TH2F("hPedHiRMS" ,"R.M.S. of pedestals, high gain",
65	rowMax,0.,rowMax, colMax,0.,colMax);
66	TH2F *hPedHiNum = new TH2F("hPedHiNum" ,"Number of pedestals, high gain",
67	rowMax,0.,rowMax, colMax,0.,colMax);
68	TH2F *hPedLoMean = new TH2F("hPedLoMean","Mean pedestals, low gain",
69	rowMax,0.,rowMax, colMax,0.,colMax);
70	TH2F *hPedLoRMS = new TH2F("hPedLoRMS" ,"R.M.S. of pedestals, low gain",
71	rowMax,0.,rowMax, colMax,0.,colMax);
72	TH2F *hPedLoNum = new TH2F("hPedLoNum" ,"Number of pedestals, low gain",
73	rowMax,0.,rowMax, colMax,0.,colMax);
74
75	TH2F *hLedHiMean = new TH2F("hLedHiMean","Mean LED, high gain",
76	rowMax,0.,rowMax, colMax,0.,colMax);
77	TH2F *hLedHiRMS = new TH2F("hLedHiRMS" ,"R.M.S. of LED, high gain",
78	rowMax,0.,rowMax, colMax,0.,colMax);
79	TH2F *hLedHiNum = new TH2F("hLedHiNum" ,"Number of LED, high gain",
80	rowMax,0.,rowMax, colMax,0.,colMax);
81	TH2F *hLedLoMean = new TH2F("hLedLoMean","Mean LED, low gain",
82	rowMax,0.,rowMax, colMax,0.,colMax);
83	TH2F *hLedLoRMS = new TH2F("hLedLoRMS" ,"R.M.S. of LED, low gain",
84	rowMax,0.,rowMax, colMax,0.,colMax);
85	TH2F *hLedLoNum = new TH2F("hLedLoNum" ,"Number of LED, low gain",
86	rowMax,0.,rowMax, colMax,0.,colMax);
87
88	Int_t iEvent=0;
89	Int_t iBin=0;
90	Int_t gain=-1,mod=-1,col=-1,row=-1;
91	Int_t runNum = 0;
92	Double_t maxAmp = 0;
93	// Double_t meanPed,rmsPed;
94	Float_t pedMean = 0;
95	Int_t nPed = 0;
96	while (rf->NextEvent()) {
97	cout << "\n\n\tEvent "<< iEvent << " of type " << rf->GetType() << endl;
98	runNum = rf->GetRunNumber();
99	while ( in.Next() ) {
100	if (!in.IsNewHWAddress() && debug > 1) {
101	cout << " time=" << in.GetTime() << "/"<< in.GetTimeLength()
102	<< " amp=" << in.GetSignal() <<endl;
103	}
104	if (in.IsNewHWAddress()) {
105	iBin = 0;
106	if (gain!=-1 && col!=-1 && row!=-1) {
107	if (nPed < 1) continue;
108	maxAmp -= (Double_t)(pedMean/nPed); // "pedestal subtraction"
109	if (maxAmp > 5) hLed[gain][col][row]->Fill(maxAmp);
110	}
111	maxAmp = 0;
112	pedMean = 0;
113	nPed = 0;
114	gain = in.IsLowGain(); // 0-high, 1-low
115	mod = in.GetModule();
116	col = in.GetColumn();
117	row = in.GetRow();
118	if (debug > 0) {
119	cout << "\n\tNew HW address: "<<in.GetHWAddress()
120	<< ", gain="<< gain
121	<< ", module="<<mod
122	<< ", xz=("<<col<<","<<row<<")"<< endl;
123	}
124
125	// Create a new histo for the new cell, if it does not exist yet
126	if (!hPed[gain][col][row]) {
127	TString name = baseNamePed;
128	TString title = baseTitlePed;
129	name +="_g"; name +=gain;
130	name +="_m"; name +=mod;
131	name +="_z"; name +=col;
132	name +="_x"; name +=row;
133
134	title +=mod; title +=",";
135	title +=col; title +=",";
136	title +=row; title +="), ";
137	title +=sgain[gain]; title +=" gain";
138
139	hPed[gain][col][row] = new TH1F(name,title,100,0.,100.);
140	hPed[gain][col][row]->Sumw2();
141	hPed[gain][col][row]->SetMarkerStyle(20);
142	hPed[gain][col][row]->SetOption("eph");
143	}
144	if (!hLed[gain][col][row]) {
145	TString name = baseNameLed;
146	TString title = baseTitleLed;
147	name +="_g"; name +=gain;
148	name +="_m"; name +=mod;
149	name +="_z"; name +=col;
150	name +="_x"; name +=row;
151
152	title +=mod; title +=",";
153	title +=col; title +=",";
154	title +=row; title +="), ";
155	title +=sgain[gain]; title +=" gain";
156
157	hLed[gain][col][row] = new TH1F(name,title,1023,0.,1023.);
158	hLed[gain][col][row]->Sumw2();
159	hLed[gain][col][row]->SetMarkerStyle(20);
160	hLed[gain][col][row]->SetOption("eph");
161	}
162	}
163
164	iBin++;
165	// Assume that first nPresample are pedestals
166	if (in.GetTimeLength()-iBin < nPresample) {
167	hPed[gain][col][row]->Fill(in.GetSignal());
168	pedMean += in.GetSignal();
169	nPed++;
170	}
171	if (in.GetSignal() > maxAmp) maxAmp = in.GetSignal();
172	} // end of next signal
173	iEvent++;
174	} // end of next event
175
176	// Fill 2-dim histograms for mean, rms and n pedestals
177
178	gain = 0;
179	for (Int_t col=0; col<colMax; col++) {
180	for (Int_t row=0; row<rowMax; row++) {
181	if (hPed[gain][col][row] != 0) {
182	hPedHiMean1->Fill( hPed[gain][col][row]->GetMean());
183	hPedHiRMS1 ->Fill( hPed[gain][col][row]->GetRMS() );
184	hPedHiMean ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetMean() );
185	hPedHiRMS ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetRMS() );
186	hPedHiNum ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetEntries() );
187	}
188	if (hLed[gain][col][row] != 0) {
189	hLedHiMean ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetMean() );
190	hLedHiRMS ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetRMS() );
191	hLedHiNum ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetEntries() );
192	}
193	}
194	}
195	gain = 1;
196	for (Int_t col=0; col<colMax; col++) {
197	for (Int_t row=0; row<rowMax; row++) {
198	if (hPed[gain][col][row] != 0) {
199	hPedLoMean1->Fill( hPed[gain][col][row]->GetMean());
200	hPedLoRMS1 ->Fill( hPed[gain][col][row]->GetRMS() );
201	hPedLoMean ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetMean() );
202	hPedLoRMS ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetRMS() );
203	hPedLoNum ->Fill( (Double_t)row, (Double_t)col, hPed[gain][col][row]->GetEntries() );
204	}
205	if (hLed[gain][col][row] != 0) {
206	hLedLoMean ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetMean() );
207	hLedLoRMS ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetRMS() );
208	hLedLoNum ->Fill( (Double_t)row, (Double_t)col, hLed[gain][col][row]->GetEntries() );
209	}
210	}
211	}
212
213	// Write existing histograms to a root file
214
215	TString fileName = "ped";
216	fileName += runNum;
217	fileName += ".root";
218	TFile *file = new TFile(fileName,"RECREATE");
219
220	for (Int_t gain=0; gain<gainMax; gain++) {
221	for (Int_t col=0; col<colMax; col++) {
222	for (Int_t row=0; row<rowMax; row++) {
223	if (hPed[gain][col][row] != 0)
224	hPed[gain][col][row]->Write();
225	if (hLed[gain][col][row] != 0)
226	hLed[gain][col][row]->Write();
227	}
228	}
229	}
230	hPedHiMean1->Write();
231	hPedHiRMS1 ->Write();
232	hPedLoMean1->Write();
233	hPedLoRMS1 ->Write();
234	hPedHiMean ->Write();
235	hPedHiRMS ->Write();
236	hPedHiNum ->Write();
237	hPedLoMean ->Write();
238	hPedLoRMS ->Write();
239	hPedLoNum ->Write();
240
241	hLedHiMean ->Write();
242	hLedHiRMS ->Write();
243	hLedHiNum ->Write();
244	hLedLoMean ->Write();
245	hLedLoRMS ->Write();
246	hLedLoNum ->Write();
247
248	file->Close();
249	stopwatch.Print();
250	}