45549361afff7140ad61d9f06454b0572f9e16a0
[u/mrichter/AliRoot.git] / TPC / AliTPCdataQA.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 /* $Id$ */
18
19
20 // stl includes
21 #include <iostream>
22
23 using namespace std;
24
25 //Root includes
26 #include <TH1F.h>
27 #include <TH2F.h>
28 #include <TString.h>
29 #include <TMath.h>
30 #include <TF1.h>
31 #include <TRandom.h>
32 #include <TDirectory.h>
33 #include <TFile.h>
34 //AliRoot includes
35 #include "AliRawReader.h"
36 #include "AliRawReaderRoot.h"
37 #include "AliRawReaderDate.h"
38 #include "AliTPCRawStream.h"
39 #include "AliTPCCalROC.h"
40 #include "AliTPCROC.h"
41 #include "AliMathBase.h"
42 #include "TTreeStream.h"
43 #include "AliTPCRawStreamFast.h"
44
45 //date
46 #include "event.h"
47 #include "AliTPCCalPad.h"
48
49 //header file
50 #include "AliTPCdataQA.h"
51
52 ClassImp(AliTPCdataQA)
53
54 AliTPCdataQA::AliTPCdataQA() : /*FOLD00*/
55   TObject(),
56   fFirstTimeBin(60),
57   fLastTimeBin(1000),
58   fAdcMin(1),
59   fAdcMax(100),
60   fOldRCUformat(kTRUE),
61   fROC(AliTPCROC::Instance()),
62   fMapping(NULL),
63   fPedestal(0),
64   fNoise(0),
65   fMaxCharge(0),
66   fMeanCharge(0),
67   fNoThreshold(0),
68   fOverThreshold0(0),
69   fOverThreshold5(0),
70   fOverThreshold10(0),
71   fOverThreshold20(0),
72   fOverThreshold30(0),
73   fEventCounter(0)
74 {
75   //
76   // default constructor
77   //
78
79   fSectorLast  = -1;
80   fRowLast     =  0;
81   fPadLast     =  0;
82   fTimeBinLast =  0;
83   fSignalLast  =  0;
84   fNAboveThreshold = 0;
85 }
86
87
88 //_____________________________________________________________________
89 AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/
90   TObject(ped),
91   fFirstTimeBin(ped.GetFirstTimeBin()),
92   fLastTimeBin(ped.GetLastTimeBin()),
93   fAdcMin(ped.GetAdcMin()),
94   fAdcMax(ped.GetAdcMax()),
95   fOldRCUformat(ped.fOldRCUformat),
96   fROC(AliTPCROC::Instance()),
97   fMapping(NULL)
98 {
99   //
100   // copy constructor
101   //
102  
103 }
104
105
106 //_____________________________________________________________________
107 AliTPCdataQA& AliTPCdataQA::operator = (const  AliTPCdataQA &source)
108 {
109   //
110   // assignment operator
111   //
112   if (&source == this) return *this;
113   new (this) AliTPCdataQA(source);
114
115   return *this;
116 }
117
118
119 //_____________________________________________________________________
120 AliTPCdataQA::~AliTPCdataQA() /*FOLD00*/
121 {
122   //
123   // destructor
124   //
125
126   // do not delete fMapping, because we do not own it.
127
128 }
129
130
131
132
133 //_____________________________________________________________________
134 Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
135 {
136   //
137   // Event Processing loop - AliTPCRawStream
138   //
139   Bool_t withInput = kFALSE;
140
141   while ( rawStreamFast->NextDDL() ){
142       while ( rawStreamFast->NextChannel() ){
143           Int_t isector  = rawStreamFast->GetSector();                       //  current sector
144           Int_t iRow     = rawStreamFast->GetRow();                          //  current row
145           Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
146
147           while ( rawStreamFast->NextBunch() ){
148             Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
149             Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
150
151               for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
152                   Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
153                   Update(isector,iRow,iPad,iTimeBin+1,signal);
154                   withInput = kTRUE;
155               }
156           }
157       }
158   }
159
160   return withInput;
161 }
162 //_____________________________________________________________________
163 Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *rawReader)
164 {
165   //
166   //  Event processing loop - AliRawReader
167   //
168   fSectorLast  = -1;
169   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
170   Bool_t res=ProcessEventFast(rawStreamFast);
171   if(res)
172     fEventCounter++; // only increment event counter if there is TPC data
173                      // otherwise Analyse (called in QA) fails
174
175   delete rawStreamFast;
176   return res;
177 }
178
179 //_____________________________________________________________________
180 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream)
181 {
182   //
183   // Event Processing loop - AliTPCRawStream
184   //
185
186   rawStream->SetOldRCUFormat(fOldRCUformat);
187
188   Bool_t withInput = kFALSE;
189
190   while (rawStream->Next()) {
191
192     Int_t iSector  = rawStream->GetSector();      //  current ROC
193     Int_t iRow     = rawStream->GetRow();         //  current row
194     Int_t iPad     = rawStream->GetPad();         //  current pad
195     Int_t iTimeBin = rawStream->GetTime();        //  current time bin
196     Float_t signal = rawStream->GetSignal();      //  current ADC signal
197     
198     Update(iSector,iRow,iPad,iTimeBin,signal);
199     withInput = kTRUE;
200   }
201
202   return withInput;
203 }
204
205
206 //_____________________________________________________________________
207 Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader)
208 {
209   //
210   //  Event processing loop - AliRawReader
211   //
212
213   // if fMapping is NULL the rawstream will crate its own mapping
214   fSectorLast  = -1;
215   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
216   rawReader->Select("TPC");
217   Bool_t res =  ProcessEvent(&rawStream);
218
219   if(res)
220     fEventCounter++; // only increment event counter if there is TPC data
221                      // otherwise Analyse (called in QA) fails
222   return res;
223 }
224
225
226 //_____________________________________________________________________
227 Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *event)
228 {
229   //
230   //  process date event
231   //
232
233   AliRawReader *rawReader = new AliRawReaderDate((void*)event);
234   Bool_t result=ProcessEvent(rawReader);
235   delete rawReader;
236   return result;
237 }
238
239
240
241 //_____________________________________________________________________
242 void AliTPCdataQA::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
243 {
244   //
245   //  Write class to file
246   //
247
248   TString sDir(dir);
249   TString option;
250
251   if ( append )
252     option = "update";
253   else
254     option = "recreate";
255
256   TDirectory *backup = gDirectory;
257   TFile f(filename,option.Data());
258   f.cd();
259   if ( !sDir.IsNull() ){
260     f.mkdir(sDir.Data());
261     f.cd(sDir);
262   }
263   this->Write();
264   f.Close();
265
266   if ( backup ) backup->cd();
267 }
268
269
270 //_____________________________________________________________________
271 Int_t AliTPCdataQA::Update(const Int_t icsector, /*FOLD00*/
272                                   const Int_t icRow,
273                                   const Int_t icPad,
274                                   const Int_t icTimeBin,
275                                   const Float_t csignal)
276 {
277   //
278   // Signal filling method
279   //
280   if (icTimeBin<fFirstTimeBin) return 0;
281   if (icTimeBin>fLastTimeBin) return 0;
282
283   if (!fMaxCharge) fMaxCharge = new AliTPCCalPad("MaxCharge","MaxCharge");
284   if (!fMeanCharge) fMeanCharge = new AliTPCCalPad("MeanCharge","MeanCharge");
285   if (!fNoThreshold) fNoThreshold = new AliTPCCalPad("NoThreshold","NoThreshold");
286   if (!fOverThreshold0) fOverThreshold0 = new AliTPCCalPad("OverThreshold0","OverThreshold0");
287   if (!fOverThreshold5) fOverThreshold5 = new AliTPCCalPad("OverThreshold5","OverThreshold5");
288   if (!fOverThreshold10) fOverThreshold10 = new AliTPCCalPad("OverThreshold10","OverThreshold10");
289   if (!fOverThreshold20) fOverThreshold20 = new AliTPCCalPad("OverThreshold20","OverThreshold20");
290   if (!fOverThreshold30) fOverThreshold30 = new AliTPCCalPad("OverThreshold30","OverThreshold30");
291   //
292
293   Int_t signal = Int_t(csignal);
294
295   // if pedestal calibrations are loaded subtract pedestals
296   if(fPedestal) {
297
298     Int_t pedestal = Int_t(fPedestal->GetCalROC(icsector)->GetValue(icRow, icPad));
299     if(pedestal<10 || pedestal>90)
300       return 0;
301     signal -= pedestal;
302   }
303
304
305   if (signal >= 0) {
306     
307     Float_t count = fNoThreshold->GetCalROC(icsector)->GetValue(icRow, icPad);
308     fNoThreshold->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
309   }
310
311   // Require at least 3 ADC channels
312   if (signal < 3)
313     return 0;
314
315   // if noise calibrations are loaded require at least 3*sigmaNoise
316   if(fNoise) {
317   
318     Float_t noise = fNoise->GetCalROC(icsector)->GetValue(icRow, icPad);
319
320     if(signal<noise*3)
321       return 0;
322   }
323   //
324   // this signal is ok - now see if the previous signal was connected
325   // this is a bit ugly since the standard decoder goes down in time bins
326   // (10, 9, 8..) while the fast HLT decoder goes up in time bins (1, 2, 3..) 
327   //
328   if(fSectorLast==icsector && fRowLast==icRow && fPadLast==icPad &&
329      fTimeBinLast==icTimeBin+1 || fTimeBinLast==icTimeBin-1)
330     fNAboveThreshold++;
331   else
332     fNAboveThreshold = 1;
333     
334   if(fNAboveThreshold==2) {
335     
336     //
337     // This is the only special case, because before we did not know if we
338     // should store the information
339     //
340     UpdateSignalHistograms(fSectorLast, fRowLast, fPadLast, fTimeBinLast,
341                          fSignalLast);
342   }
343   
344   // keep the information for the next signal
345   fSectorLast  = icsector;
346   fRowLast     = icRow;
347   fPadLast     = icPad;
348   fTimeBinLast = icTimeBin;
349   fSignalLast  = signal;
350   
351   if(fNAboveThreshold==1) // we don't know if this information should be stored
352     return 1;
353   
354   UpdateSignalHistograms(fSectorLast, fRowLast, fPadLast, fTimeBinLast,
355                        fSignalLast);
356
357   return 1;
358 }
359 //_____________________________________________________________________
360 void AliTPCdataQA::UpdateSignalHistograms(const Int_t icsector, /*FOLD00*/
361                                         const Int_t icRow,
362                                         const Int_t icPad,
363                                         const Int_t icTimeBin,
364                                         const Float_t signal)
365 {
366   //
367   // Signal filling method
368   //
369   
370   {
371     Float_t charge = fMeanCharge->GetCalROC(icsector)->GetValue(icRow, icPad);
372     fMeanCharge->GetCalROC(icsector)->SetValue(icRow, icPad, charge + signal);
373   }
374   
375   if (signal>fMaxCharge->GetCalROC(icsector)->GetValue(icRow, icPad)){
376     fMaxCharge->GetCalROC(icsector)->SetValue(icRow, icPad,signal);
377   }
378   
379   if (signal>0){
380     Float_t count = fOverThreshold0->GetCalROC(icsector)->GetValue(icRow, icPad);
381     fOverThreshold0->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
382   };
383   //
384   if (signal>5){
385     Float_t count = fOverThreshold5->GetCalROC(icsector)->GetValue(icRow, icPad);
386     fOverThreshold5->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
387   };
388   if (signal>10){
389     Float_t count = fOverThreshold10->GetCalROC(icsector)->GetValue(icRow, icPad);
390     fOverThreshold10->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
391   };
392   if (signal>20){
393     Float_t count = fOverThreshold20->GetCalROC(icsector)->GetValue(icRow, icPad);
394     fOverThreshold20->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
395   };
396   if (signal>30){
397     Float_t count = fOverThreshold30->GetCalROC(icsector)->GetValue(icRow, icPad);
398     fOverThreshold30->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
399   };  
400 }
401
402 //_____________________________________________________________________
403 void AliTPCdataQA::Analyse()
404 {
405   //
406   //  Calculate calibration constants
407   //
408   
409   cout << "Analyse called" << endl;
410
411   if(fEventCounter==0) {
412
413     cout << "EventCounter == 0, Cannot analyse" << endl;
414     return;
415   }
416
417   Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
418   
419   cout << "EventCounter: " << fEventCounter << endl
420        << "TimeBins: " << nTimeBins << endl;
421
422   if (fMeanCharge && fNoThreshold) fMeanCharge->Divide(fNoThreshold);
423
424   Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins);
425   if (fNoThreshold)     fNoThreshold->Multiply(normalization);  
426   if (fOverThreshold0)  fOverThreshold0->Multiply(normalization);  
427   if (fOverThreshold5)  fOverThreshold5->Multiply(normalization);  
428   if (fOverThreshold10) fOverThreshold10->Multiply(normalization);  
429   if (fOverThreshold20) fOverThreshold20->Multiply(normalization);  
430   if (fOverThreshold30) fOverThreshold30->Multiply(normalization);  
431 }