2e9ff87f41c007716ab36dc27c6047b08e21401b
[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 #include "AliTPCPreprocessorOnline.h"
49
50 //header file
51 #include "AliTPCdataQA.h"
52
53 ClassImp(AliTPCdataQA)
54
55 AliTPCdataQA::AliTPCdataQA() : /*FOLD00*/  
56   TH1F("TPCRAW","TPCRAW",100,0,100),
57   fFirstTimeBin(60),
58   fLastTimeBin(1000),
59   fMaxTime(1100),
60   fAdcMin(1),
61   fAdcMax(100),
62   fOldRCUformat(kTRUE),
63   fROC(AliTPCROC::Instance()),
64   fMapping(NULL),
65   fPedestal(0),
66   fNoise(0),
67   fMaxCharge(0),
68   fMeanCharge(0),
69   fNoThreshold(0),
70   fOverThreshold0(0),
71   fOverThreshold5(0),
72   fOverThreshold10(0),
73   fOverThreshold20(0),
74   fOverThreshold30(0),
75   fEventCounter(0)
76 {
77   //
78   // default constructor
79   //
80
81   fSectorLast  = -1;
82   fRowLast     =  0;
83   fPadLast     =  0;
84   fTimeBinLast =  0;
85   fSignalLast  =  0;
86   fNAboveThreshold = 0;
87 }
88
89
90 //_____________________________________________________________________
91 AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/
92   TH1F(ped),
93   fFirstTimeBin(ped.GetFirstTimeBin()),
94   fLastTimeBin(ped.GetLastTimeBin()),
95   fMaxTime(ped.fMaxTime),
96   fAdcMin(ped.GetAdcMin()),
97   fAdcMax(ped.GetAdcMax()),
98   fOldRCUformat(ped.fOldRCUformat),
99   fROC(AliTPCROC::Instance()),
100   fMapping(NULL)
101 {
102   //
103   // copy constructor
104   //
105  
106 }
107
108
109 //_____________________________________________________________________
110 AliTPCdataQA& AliTPCdataQA::operator = (const  AliTPCdataQA &source)
111 {
112   //
113   // assignment operator
114   //
115   if (&source == this) return *this;
116   new (this) AliTPCdataQA(source);
117
118   return *this;
119 }
120
121
122 //_____________________________________________________________________
123 AliTPCdataQA::~AliTPCdataQA() /*FOLD00*/
124 {
125   //
126   // destructor
127   //
128
129   // do not delete fMapping, because we do not own it.
130
131 }
132
133
134
135
136 //_____________________________________________________________________
137 Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
138 {
139   //
140   // Event Processing loop - AliTPCRawStream
141   //
142   Bool_t withInput = kFALSE;
143
144   while ( rawStreamFast->NextDDL() ){
145       while ( rawStreamFast->NextChannel() ){
146           Int_t isector  = rawStreamFast->GetSector();                       //  current sector
147           Int_t iRow     = rawStreamFast->GetRow();                          //  current row
148           Int_t iPad     = rawStreamFast->GetPad();                          //  current pad
149
150           while ( rawStreamFast->NextBunch() ){
151             Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
152             Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
153
154               for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
155                   Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
156                   Update(isector,iRow,iPad,iTimeBin+1,signal);
157                   withInput = kTRUE;
158               }
159           }
160       }
161   }
162
163   return withInput;
164 }
165 //_____________________________________________________________________
166 Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *rawReader)
167 {
168   //
169   //  Event processing loop - AliRawReader
170   //
171   fSectorLast  = -1;
172   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
173   Bool_t res=ProcessEventFast(rawStreamFast);
174   if(res)
175     fEventCounter++; // only increment event counter if there is TPC data
176                      // otherwise Analyse (called in QA) fails
177
178   delete rawStreamFast;
179   return res;
180 }
181
182 //_____________________________________________________________________
183 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream)
184 {
185   //
186   // Event Processing loop - AliTPCRawStream
187   //
188
189   rawStream->SetOldRCUFormat(fOldRCUformat);
190
191   Bool_t withInput = kFALSE;
192
193   while (rawStream->Next()) {
194
195     Int_t iSector  = rawStream->GetSector();      //  current ROC
196     Int_t iRow     = rawStream->GetRow();         //  current row
197     Int_t iPad     = rawStream->GetPad();         //  current pad
198     Int_t iTimeBin = rawStream->GetTime();        //  current time bin
199     Float_t signal = rawStream->GetSignal();      //  current ADC signal
200     
201     Update(iSector,iRow,iPad,iTimeBin,signal);
202     withInput = kTRUE;
203   }
204
205   return withInput;
206 }
207
208
209 //_____________________________________________________________________
210 Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader)
211 {
212   //
213   //  Event processing loop - AliRawReader
214   //
215
216   // if fMapping is NULL the rawstream will crate its own mapping
217   fSectorLast  = -1;
218   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
219   rawReader->Select("TPC");
220   Bool_t res =  ProcessEvent(&rawStream);
221
222   if(res)
223     fEventCounter++; // only increment event counter if there is TPC data
224                      // otherwise Analyse (called in QA) fails
225   return res;
226 }
227
228
229 //_____________________________________________________________________
230 Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *event)
231 {
232   //
233   //  process date event
234   //
235
236   AliRawReader *rawReader = new AliRawReaderDate((void*)event);
237   Bool_t result=ProcessEvent(rawReader);
238   delete rawReader;
239   return result;
240 }
241
242
243
244 //_____________________________________________________________________
245 void AliTPCdataQA::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
246 {
247   //
248   //  Write class to file
249   //
250
251   TString sDir(dir);
252   TString option;
253
254   if ( append )
255     option = "update";
256   else
257     option = "recreate";
258
259   TDirectory *backup = gDirectory;
260   TFile f(filename,option.Data());
261   f.cd();
262   if ( !sDir.IsNull() ){
263     f.mkdir(sDir.Data());
264     f.cd(sDir);
265   }
266   this->Write();
267   f.Close();
268
269   if ( backup ) backup->cd();
270 }
271
272
273 //_____________________________________________________________________
274 Int_t AliTPCdataQA::Update(const Int_t icsector, /*FOLD00*/
275                                   const Int_t icRow,
276                                   const Int_t icPad,
277                                   const Int_t icTimeBin,
278                                   const Float_t csignal)
279 {
280   //
281   // Signal filling method
282   //
283   if (icTimeBin<fFirstTimeBin) return 0;
284   if (icTimeBin>fLastTimeBin) return 0;
285
286   if (!fMaxCharge) fMaxCharge = new AliTPCCalPad("MaxCharge","MaxCharge");
287   if (!fMeanCharge) fMeanCharge = new AliTPCCalPad("MeanCharge","MeanCharge");
288   if (!fNoThreshold) fNoThreshold = new AliTPCCalPad("NoThreshold","NoThreshold");
289   if (!fOverThreshold0) fOverThreshold0 = new AliTPCCalPad("OverThreshold0","OverThreshold0");
290   if (!fOverThreshold5) fOverThreshold5 = new AliTPCCalPad("OverThreshold5","OverThreshold5");
291   if (!fOverThreshold10) fOverThreshold10 = new AliTPCCalPad("OverThreshold10","OverThreshold10");
292   if (!fOverThreshold20) fOverThreshold20 = new AliTPCCalPad("OverThreshold20","OverThreshold20");
293   if (!fOverThreshold30) fOverThreshold30 = new AliTPCCalPad("OverThreshold30","OverThreshold30");
294   //
295
296   Int_t signal = Int_t(csignal);
297
298   // if pedestal calibrations are loaded subtract pedestals
299   if(fPedestal) {
300
301     Int_t pedestal = Int_t(fPedestal->GetCalROC(icsector)->GetValue(icRow, icPad));
302     if(pedestal<10 || pedestal>90)
303       return 0;
304     signal -= pedestal;
305   }
306
307
308   if (signal >= 0) {
309     
310     Float_t count = fNoThreshold->GetCalROC(icsector)->GetValue(icRow, icPad);
311     fNoThreshold->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
312   }
313
314   // Require at least 3 ADC channels
315   if (signal < 3)
316     return 0;
317
318   // if noise calibrations are loaded require at least 3*sigmaNoise
319   if(fNoise) {
320   
321     Float_t noise = fNoise->GetCalROC(icsector)->GetValue(icRow, icPad);
322
323     if(signal<noise*3)
324       return 0;
325   }
326   //
327   // this signal is ok - now see if the previous signal was connected
328   // this is a bit ugly since the standard decoder goes down in time bins
329   // (10, 9, 8..) while the fast HLT decoder goes up in time bins (1, 2, 3..) 
330   //
331   if(fSectorLast==icsector && fRowLast==icRow && fPadLast==icPad &&
332      fTimeBinLast==icTimeBin+1 || fTimeBinLast==icTimeBin-1)
333     fNAboveThreshold++;
334   else
335     fNAboveThreshold = 1;
336     
337   if(fNAboveThreshold==2) {
338     
339     //
340     // This is the only special case, because before we did not know if we
341     // should store the information
342     //
343     UpdateSignalHistograms(fSectorLast, fRowLast, fPadLast, fTimeBinLast,
344                          fSignalLast);
345   }
346   
347   // keep the information for the next signal
348   fSectorLast  = icsector;
349   fRowLast     = icRow;
350   fPadLast     = icPad;
351   fTimeBinLast = icTimeBin;
352   fSignalLast  = signal;
353   
354   if(fNAboveThreshold==1) // we don't know if this information should be stored
355     return 1;
356   
357   UpdateSignalHistograms(fSectorLast, fRowLast, fPadLast, fTimeBinLast,
358                        fSignalLast);
359
360   return 1;
361 }
362 //_____________________________________________________________________
363 void AliTPCdataQA::UpdateSignalHistograms(const Int_t icsector, /*FOLD00*/
364                                         const Int_t icRow,
365                                         const Int_t icPad,
366                                         const Int_t icTimeBin,
367                                         const Float_t signal)
368 {
369   //
370   // Signal filling method
371   //
372   
373   {
374     Float_t charge = fMeanCharge->GetCalROC(icsector)->GetValue(icRow, icPad);
375     fMeanCharge->GetCalROC(icsector)->SetValue(icRow, icPad, charge + signal);
376   }
377   
378   if (signal>fMaxCharge->GetCalROC(icsector)->GetValue(icRow, icPad)){
379     fMaxCharge->GetCalROC(icsector)->SetValue(icRow, icPad,signal);
380   }
381   
382   if (signal>0){
383     Float_t count = fOverThreshold0->GetCalROC(icsector)->GetValue(icRow, icPad);
384     fOverThreshold0->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
385   };
386   //
387   if (signal>5){
388     Float_t count = fOverThreshold5->GetCalROC(icsector)->GetValue(icRow, icPad);
389     fOverThreshold5->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
390   };
391   if (signal>10){
392     Float_t count = fOverThreshold10->GetCalROC(icsector)->GetValue(icRow, icPad);
393     fOverThreshold10->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
394   };
395   if (signal>20){
396     Float_t count = fOverThreshold20->GetCalROC(icsector)->GetValue(icRow, icPad);
397     fOverThreshold20->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
398   };
399   if (signal>30){
400     Float_t count = fOverThreshold30->GetCalROC(icsector)->GetValue(icRow, icPad);
401     fOverThreshold30->GetCalROC(icsector)->SetValue(icRow, icPad,count+1);
402   };  
403 }
404
405 //_____________________________________________________________________
406 void AliTPCdataQA::Analyse()
407 {
408   //
409   //  Calculate calibration constants
410   //
411   
412   cout << "Analyse called" << endl;
413
414   if(fEventCounter==0) {
415
416     cout << "EventCounter == 0, Cannot analyse" << endl;
417     return;
418   }
419
420   Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
421   cout << "EventCounter: " << fEventCounter << endl
422        << "TimeBins: " << nTimeBins << endl;
423
424   if (fMeanCharge && fNoThreshold) fMeanCharge->Divide(fNoThreshold);
425
426   Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins);
427   if (fNoThreshold)     fNoThreshold->Multiply(normalization);  
428   if (fOverThreshold0)  fOverThreshold0->Multiply(normalization);  
429   if (fOverThreshold5)  fOverThreshold5->Multiply(normalization);  
430   if (fOverThreshold10) fOverThreshold10->Multiply(normalization);  
431   if (fOverThreshold20) fOverThreshold20->Multiply(normalization);  
432   if (fOverThreshold30) fOverThreshold30->Multiply(normalization);  
433 }
434
435
436 void AliTPCdataQA::MakeTree(const char *fname){
437   //
438   // Export result to the tree -located in the file
439   // This file can be analyzed using AliTPCCalibViewer
440   // 
441   AliTPCdataQA *ped = this;
442   AliTPCPreprocessorOnline preprocesor;
443   if (ped->GetMaxCharge()) preprocesor.AddComponent(ped->GetMaxCharge());  
444   if (ped->GetMeanCharge()) preprocesor.AddComponent(ped->GetMeanCharge());  
445   if (ped->GetOverThreshold0()) preprocesor.AddComponent(ped->GetOverThreshold0());
446   if (ped->GetOverThreshold5()) preprocesor.AddComponent(ped->GetOverThreshold5());
447   if (ped->GetOverThreshold10()) preprocesor.AddComponent(ped->GetOverThreshold10());
448   if (ped->GetOverThreshold20()) preprocesor.AddComponent(ped->GetOverThreshold20());
449   if (ped->GetOverThreshold30()) preprocesor.AddComponent(ped->GetOverThreshold30());
450   preprocesor.DumpToFile(fname);  
451 }
452
453
454
455 void AliTPCdataQA::MakeArrays(){
456   //
457   //
458   //
459   AliTPCROC * roc = AliTPCROC::Instance();
460   //
461   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
462   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1);
463  
464   fAllBins = new Float_t*[nRowsMax];
465   fAllSigBins = new Int_t*[nRowsMax];
466   fAllNSigBins = new Int_t[nRowsMax];
467   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
468     //
469     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
470     fAllBins[iRow] = new Float_t[maxBin];
471     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
472     fAllSigBins[iRow] = new Int_t[maxBin];
473     fAllNSigBins[iRow]=0;
474   }
475 }
476
477
478 void AliTPCdataQA::CleanArrays(){
479   //
480   //
481   //
482   AliTPCROC * roc = AliTPCROC::Instance();
483   //
484   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
485   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1); 
486   for (Int_t iRow = 0; iRow < nRowsMax; iRow++) {
487     //
488     Int_t maxBin = fMaxTime*(nPadsMax+6);  // add 3 virtual pads  before and 3 after
489     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
490     fAllNSigBins[iRow]=0;
491   }
492 }
493
494 Float_t* AliTPCdataQA::GetExpandDigit(Int_t row, Int_t pad, Int_t time){
495   //
496   //
497   //
498   AliTPCROC * roc = AliTPCROC::Instance();
499   Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1);
500   if (row<0 || row>nRowsMax) return 0;
501   Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1); 
502   
503 }