Example of low-level data processing
[u/mrichter/AliRoot.git] / PHOS / macros / RawProcessing / PedAmp.C
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 }