bfaf56c66a52400b7d10634db6ebae6b2ab7a8e5
[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   February 2008
21
22   The code has been heavily modified so that now the RAW data is
23   "expanded" for each sector and stored in a big signal array. Then a
24   simple version of the code in AliTPCclustererMI is used to identify
25   the local maxima and these are then used for QA. This gives a better
26   estimate of the charge (both max and total) and also limits the
27   effect of noise.
28
29   Implementation:
30
31   In Update the RAW signals >= 3 ADC channels are stored in the arrays.
32   
33   There are 3 arrays:
34   Float_t** fAllBins       2d array [row][bin(pad, time)] ADC signal
35   Int_t**   fAllSigBins    2d array [row][signal#] bin(with signal)
36   Int_t*    fAllNSigBins;  1d array [row] Nsignals
37
38   This is done sector by sector.
39
40   When data from a new sector is encountered, the method
41   FindLocalMaxima is called on the data from the previous sector, and
42   the calibration/data objects are updated with the "cluster"
43   info. Finally the arrays are cleared.
44
45   The requirements for a local maxima is:
46   Charge in bin is >= 5 ADC channels.
47   Charge in bin is larger than all the 8 neighboring bins.
48   At least one of the two pad neighbors has a signal.
49   At least one of the two time neighbors has a signal.
50
51   Before accessing the data it is expected that the Analyse method is
52   called. This normalizes some of the data objects to per event or per
53   cluster. 
54   If more data is passed to the class after Analyse has been called
55   the normalization is reversed and Analyse has to be called again.
56 */
57
58
59 // stl includes
60 #include <iostream>
61
62 using namespace std;
63
64 //Root includes
65 #include <TH1F.h>
66 #include <TH2F.h>
67 #include <TString.h>
68 #include <TMath.h>
69 #include <TF1.h>
70 #include <TRandom.h>
71 #include <TDirectory.h>
72 #include <TFile.h>
73 #include <TError.h>
74 #include <TMap.h>
75 //AliRoot includes
76 #include "AliRawReader.h"
77 #include "AliRawReaderRoot.h"
78 #include "AliRawReaderDate.h"
79 #include "AliTPCRawStream.h"
80 #include "AliTPCRawStreamV3.h"
81 #include "AliTPCCalROC.h"
82 #include "AliTPCROC.h"
83 #include "AliMathBase.h"
84 #include "TTreeStream.h"
85 #include "AliTPCRawStreamFast.h"
86
87 //date
88 #include "event.h"
89 #include "AliTPCCalPad.h"
90 #include "AliTPCPreprocessorOnline.h"
91
92 //header file
93 #include "AliTPCdataQA.h"
94
95 ClassImp(AliTPCdataQA)
96
97 AliTPCdataQA::AliTPCdataQA() : /*FOLD00*/  
98   fFirstTimeBin(60),
99   fLastTimeBin(1000),
100   fAdcMin(1),
101   fAdcMax(100),
102   fMapping(NULL),
103   fPedestal(0),
104   fNoise(0),
105   fNLocalMaxima(0),
106   fMaxCharge(0),
107   fMeanCharge(0),
108   fNoThreshold(0),
109   fNTimeBins(0),
110   fNPads(0),
111   fTimePosition(0),
112   fOverThreshold10(0),
113   fOverThreshold20(0),
114   fOverThreshold30(0),
115   fHistQVsTimeSideA(0),
116   fHistQVsTimeSideC(0),
117   fHistQMaxVsTimeSideA(0),
118   fHistQMaxVsTimeSideC(0),
119   fEventCounter(0),
120   fIsAnalysed(kFALSE),
121   fAllBins(0),
122   fAllSigBins(0),
123   fAllNSigBins(0),
124   fRowsMax(0),
125   fPadsMax(0),
126   fTimeBinsMax(0)
127 {
128   //
129   // default constructor
130   //
131 }
132
133 //_____________________________________________________________________
134 AliTPCdataQA::AliTPCdataQA(AliRecoParam::EventSpecie_t es) :
135 fFirstTimeBin(60),
136 fLastTimeBin(1000),
137 fAdcMin(1),
138 fAdcMax(100),
139 fMapping(NULL),
140 fPedestal(0),
141 fNoise(0),
142 fNLocalMaxima(0),
143 fMaxCharge(0),
144 fMeanCharge(0),
145 fNoThreshold(0),
146 fNTimeBins(0),
147 fNPads(0),
148 fTimePosition(0),
149 fOverThreshold10(0),
150 fOverThreshold20(0),
151 fOverThreshold30(0),
152 fHistQVsTimeSideA(0),
153 fHistQVsTimeSideC(0),
154 fHistQMaxVsTimeSideA(0),
155 fHistQMaxVsTimeSideC(0),
156 fEventCounter(0),
157 fIsAnalysed(kFALSE),
158 fAllBins(0),
159 fAllSigBins(0),
160 fAllNSigBins(0),
161 fRowsMax(0),
162 fPadsMax(0),
163 fTimeBinsMax(0)
164 {
165 // ctor creating the histogram
166   char *  name = Form("TPCRAW_%s", AliRecoParam::GetEventSpecieName(es)) ; 
167   TH1F(name, name,100,0,100) ; 
168 }
169
170 //_____________________________________________________________________
171 AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/
172   TH1F(ped),
173   fFirstTimeBin(ped.GetFirstTimeBin()),
174   fLastTimeBin(ped.GetLastTimeBin()),
175   fAdcMin(ped.GetAdcMin()),
176   fAdcMax(ped.GetAdcMax()),
177   fMapping(NULL),
178   fPedestal(0),
179   fNoise(0),
180   fNLocalMaxima(0),
181   fMaxCharge(0),
182   fMeanCharge(0),
183   fNoThreshold(0),
184   fNTimeBins(0),
185   fNPads(0),
186   fTimePosition(0),
187   fOverThreshold10(0),
188   fOverThreshold20(0),
189   fOverThreshold30(0),
190   fHistQVsTimeSideA(0),
191   fHistQVsTimeSideC(0),
192   fHistQMaxVsTimeSideA(0),
193   fHistQMaxVsTimeSideC(0),
194   fEventCounter(ped.GetEventCounter()),
195   fIsAnalysed(ped.GetIsAnalysed()),
196   fAllBins(0),
197   fAllSigBins(0),
198   fAllNSigBins(0),
199   fRowsMax(0),
200   fPadsMax(0),
201   fTimeBinsMax(0)
202 {
203   //
204   // copy constructor
205   //
206   if(ped.GetNLocalMaxima())
207     fNLocalMaxima  = new AliTPCCalPad(*ped.GetNLocalMaxima());
208   if(ped.GetMaxCharge())
209     fMaxCharge      = new AliTPCCalPad(*ped.GetMaxCharge());
210   if(ped.GetMeanCharge())
211     fMeanCharge     = new AliTPCCalPad(*ped.GetMeanCharge());
212   if(ped.GetNoThreshold())
213     fNoThreshold  = new AliTPCCalPad(*ped.GetNoThreshold());
214   if(ped.GetNTimeBins())
215     fNTimeBins  = new AliTPCCalPad(*ped.GetNTimeBins());
216   if(ped.GetNPads())
217     fNPads  = new AliTPCCalPad(*ped.GetNPads());
218   if(ped.GetTimePosition())
219     fTimePosition  = new AliTPCCalPad(*ped.GetTimePosition());
220   if(ped.GetOverThreshold10())
221     fOverThreshold10  = new AliTPCCalPad(*ped.GetOverThreshold10());
222   if(ped.GetOverThreshold20())
223     fOverThreshold20  = new AliTPCCalPad(*ped.GetOverThreshold20());
224   if(ped.GetOverThreshold30())
225     fOverThreshold30  = new AliTPCCalPad(*ped.GetOverThreshold30());
226   if(ped.GetHistQVsTimeSideA()) {
227     fHistQVsTimeSideA = new TProfile(*ped.GetHistQVsTimeSideA());
228     fHistQVsTimeSideA->SetDirectory(0);
229   }
230   if(ped.GetHistQVsTimeSideC()) {
231     fHistQVsTimeSideC = new TProfile(*ped.GetHistQVsTimeSideC());
232     fHistQVsTimeSideC->SetDirectory(0);
233   }
234   if(ped.GetHistQMaxVsTimeSideA()) {
235     fHistQMaxVsTimeSideA = new TProfile(*ped.GetHistQMaxVsTimeSideA());
236     fHistQMaxVsTimeSideA->SetDirectory(0);
237   }
238   if(ped.GetHistQMaxVsTimeSideC()) {
239     fHistQMaxVsTimeSideC = new TProfile(*ped.GetHistQMaxVsTimeSideC());
240     fHistQMaxVsTimeSideC->SetDirectory(0);
241   }
242 }
243
244 //_____________________________________________________________________
245 AliTPCdataQA::AliTPCdataQA(const TMap *config) : /*FOLD00*/  
246   TH1F("TPCRAW","TPCRAW",100,0,100),
247   fFirstTimeBin(60),
248   fLastTimeBin(1000),
249   fAdcMin(1),
250   fAdcMax(100),
251   fMapping(NULL),
252   fPedestal(0),
253   fNoise(0),
254   fNLocalMaxima(0),
255   fMaxCharge(0),
256   fMeanCharge(0),
257   fNoThreshold(0),
258   fNTimeBins(0),
259   fNPads(0),
260   fTimePosition(0),
261   fOverThreshold10(0),
262   fOverThreshold20(0),
263   fOverThreshold30(0),
264   fHistQVsTimeSideA(0),
265   fHistQVsTimeSideC(0),
266   fHistQMaxVsTimeSideA(0),
267   fHistQMaxVsTimeSideC(0),
268   fEventCounter(0),
269   fIsAnalysed(kFALSE),
270   fAllBins(0),
271   fAllSigBins(0),
272   fAllNSigBins(0),
273   fRowsMax(0),
274   fPadsMax(0),
275   fTimeBinsMax(0)
276 {
277   //
278   // default constructor
279   //
280   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
281   if (config->GetValue("LastTimeBin"))  fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
282   if (config->GetValue("AdcMin"))       fAdcMin = ((TObjString*)config->GetValue("AdcMin"))->GetString().Atoi();
283   if (config->GetValue("AdcMax"))       fAdcMax = ((TObjString*)config->GetValue("AdcMax"))->GetString().Atoi();
284 }
285
286 //_____________________________________________________________________
287 AliTPCdataQA& AliTPCdataQA::operator = (const  AliTPCdataQA &source)
288 {
289   //
290   // assignment operator
291   //
292   if (&source == this) return *this;
293   new (this) AliTPCdataQA(source);
294
295   return *this;
296 }
297
298
299 //_____________________________________________________________________
300 AliTPCdataQA::~AliTPCdataQA() /*FOLD00*/
301 {
302   //
303   // destructor
304   //
305
306   // do not delete fMapping, because we do not own it.
307   // do not delete fMapping, because we do not own it.
308   // do not delete fNoise and fPedestal, because we do not own them.
309
310   delete fNLocalMaxima;
311   delete fMaxCharge;
312   delete fMeanCharge;
313   delete fNoThreshold;
314   delete fNTimeBins;
315   delete fNPads;
316   delete fTimePosition;
317   delete fOverThreshold10;
318   delete fOverThreshold20;
319   delete fOverThreshold30;
320   delete fHistQVsTimeSideA;
321   delete fHistQVsTimeSideC;
322   delete fHistQMaxVsTimeSideA;
323   delete fHistQMaxVsTimeSideC;
324
325   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
326     delete [] fAllBins[iRow];
327     delete [] fAllSigBins[iRow];
328   }  
329   delete [] fAllBins;
330   delete [] fAllSigBins;
331   delete [] fAllNSigBins;
332 }
333 //_____________________________________________________________________
334 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStreamV3 *rawStreamV3)
335 {
336   //
337   // Event Processing loop - AliTPCRawStreamV3
338   //
339   Bool_t withInput = kFALSE;
340   Int_t nSignals = 0;
341   Int_t lastSector = -1;
342   
343   while ( rawStreamV3->NextDDL() ){
344     while ( rawStreamV3->NextChannel() ){
345       Int_t iSector = rawStreamV3->GetSector(); //  current sector
346       Int_t iRow    = rawStreamV3->GetRow();    //  current row
347       Int_t iPad    = rawStreamV3->GetPad();    //  current pad
348       if (iRow<0 || iPad<0) continue;
349       // Call local maxima finder if the data is in a new sector
350       if(iSector != lastSector) {
351         
352         if(nSignals>0)
353           FindLocalMaxima(lastSector);
354         
355         CleanArrays();
356         lastSector = iSector;
357         nSignals = 0;
358       }
359       
360       while ( rawStreamV3->NextBunch() ){
361         Int_t  startTbin    = (Int_t)rawStreamV3->GetStartTimeBin();
362         Int_t  bunchlength  = (Int_t)rawStreamV3->GetBunchLength();
363         const UShort_t *sig = rawStreamV3->GetSignals();
364         
365         for (Int_t iTimeBin = 0; iTimeBin<bunchlength; iTimeBin++){
366           Float_t signal=(Float_t)sig[iTimeBin];
367           nSignals += Update(iSector,iRow,iPad,startTbin--,signal);
368           withInput = kTRUE;
369         }
370       }
371     }
372   }
373   
374   if (lastSector>=0&&nSignals>0)
375     FindLocalMaxima(lastSector);
376   
377   return withInput;
378 }
379
380 //_____________________________________________________________________
381 Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader)
382 {
383   //
384   //  Event processing loop - AliRawReader
385   //
386   AliTPCRawStreamV3 rawStreamV3(rawReader, (AliAltroMapping**)fMapping);
387   Bool_t res=ProcessEvent(&rawStreamV3);
388   if(res)
389     fEventCounter++; // only increment event counter if there is TPC data
390   return res;
391 }
392
393 //_____________________________________________________________________
394 Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
395 {
396   //
397   // Event Processing loop - AliTPCRawStream
398   //
399   Bool_t withInput = kFALSE;
400   Int_t nSignals = 0;
401   Int_t lastSector = -1;
402
403   while ( rawStreamFast->NextDDL() ){
404     while ( rawStreamFast->NextChannel() ){
405       
406       Int_t iSector  = rawStreamFast->GetSector(); //  current sector
407       Int_t iRow     = rawStreamFast->GetRow();    //  current row
408       Int_t iPad     = rawStreamFast->GetPad();    //  current pad
409   // Call local maxima finder if the data is in a new sector
410       if(iSector != lastSector) {
411         
412         if(nSignals>0)
413           FindLocalMaxima(lastSector);
414         
415         CleanArrays();
416         lastSector = iSector;
417         nSignals = 0;
418       }
419       
420       while ( rawStreamFast->NextBunch() ){
421         Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
422         Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
423         
424         for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
425           Float_t signal = rawStreamFast->GetSignals()[iTimeBin-startTbin];
426           nSignals += Update(iSector,iRow,iPad,iTimeBin+1,signal);
427           withInput = kTRUE;
428         }
429       }
430     }
431   }
432   
433   return withInput;
434 }
435 //_____________________________________________________________________
436 Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *rawReader)
437 {
438   //
439   //  Event processing loop - AliRawReader
440   //
441   AliTPCRawStreamFast rawStreamFast(rawReader, (AliAltroMapping**)fMapping);
442   Bool_t res=ProcessEventFast(&rawStreamFast);
443   if(res)
444     fEventCounter++; // only increment event counter if there is TPC data
445                      // otherwise Analyse (called in QA) fails
446
447   return res;
448 }
449
450 //_____________________________________________________________________
451 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream)
452 {
453   //
454   // Event Processing loop - AliTPCRawStream
455   //
456
457   Bool_t withInput = kFALSE;
458   Int_t nSignals = 0;
459   Int_t lastSector = -1;
460
461   while (rawStream->Next()) {
462
463     Int_t iSector  = rawStream->GetSector();      //  current ROC
464     Int_t iRow     = rawStream->GetRow();         //  current row
465     Int_t iPad     = rawStream->GetPad();         //  current pad
466     Int_t iTimeBin = rawStream->GetTime();        //  current time bin
467     Float_t signal = rawStream->GetSignal();      //  current ADC signal
468     
469     // Call local maxima finder if the data is in a new sector
470     if(iSector != lastSector) {
471       
472       if(nSignals>0)
473         FindLocalMaxima(lastSector);
474       
475       CleanArrays();
476       lastSector = iSector;
477       nSignals = 0;
478     }
479     
480     // Sometimes iRow==-1 if there is problems to read the data
481     if(iRow>=0) {
482       nSignals += Update(iSector,iRow,iPad,iTimeBin,signal);
483       withInput = kTRUE;
484     }
485   }
486
487   if (lastSector>=0&&nSignals>0)
488     FindLocalMaxima(lastSector);
489   
490   return withInput;
491 }
492
493
494 //_____________________________________________________________________
495 Bool_t AliTPCdataQA::ProcessEventOld(AliRawReader *rawReader)
496 {
497   //
498   //  Event processing loop - AliRawReader
499   //
500
501   // if fMapping is NULL the rawstream will crate its own mapping
502   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
503   rawReader->Select("TPC");
504   Bool_t res =  ProcessEvent(&rawStream);
505
506   if(res)
507     fEventCounter++; // only increment event counter if there is TPC data
508                      // otherwise Analyse (called in QA) fails
509   return res;
510 }
511
512
513 //_____________________________________________________________________
514 Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *event)
515 {
516   //
517   //  process date event
518   //
519
520   AliRawReaderDate rawReader((void*)event);
521   Bool_t result=ProcessEvent(&rawReader);
522   return result;
523 }
524
525
526
527 //_____________________________________________________________________
528 void AliTPCdataQA::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
529 {
530   //
531   //  Write class to file
532   //
533
534   TString sDir(dir);
535   TString option;
536
537   if ( append )
538     option = "update";
539   else
540     option = "recreate";
541
542   TDirectory *backup = gDirectory;
543   TFile f(filename,option.Data());
544   f.cd();
545   if ( !sDir.IsNull() ){
546     f.mkdir(sDir.Data());
547     f.cd(sDir);
548   }
549   this->Write();
550   f.Close();
551
552   if ( backup ) backup->cd();
553 }
554
555
556 //_____________________________________________________________________
557 Int_t AliTPCdataQA::Update(const Int_t iSector, /*FOLD00*/
558                            const Int_t iRow,
559                            const Int_t iPad,
560                            const Int_t iTimeBin,
561                            Float_t signal)
562 {
563   //
564   // Signal filling method
565   //
566   
567   //
568   // Define the calibration objects the first time Update is called
569   // NB! This has to be done first even if the data is rejected by the time 
570   // cut to make sure that the objects are available in Analyse
571   //
572   if (!fNLocalMaxima) fNLocalMaxima = new AliTPCCalPad("NLocalMaxima","NLocalMaxima");
573   if (!fMaxCharge) fMaxCharge = new AliTPCCalPad("MaxCharge","MaxCharge");
574   if (!fMeanCharge) fMeanCharge = new AliTPCCalPad("MeanCharge","MeanCharge");
575   if (!fNoThreshold) fNoThreshold = new AliTPCCalPad("NoThreshold","NoThreshold");
576   if (!fNTimeBins) fNTimeBins = new AliTPCCalPad("NTimeBins","NTimeBins");
577   if (!fNPads) fNPads = new AliTPCCalPad("NPads","NPads");
578   if (!fTimePosition) fTimePosition = new AliTPCCalPad("TimePosition","TimePosition");
579   if (!fOverThreshold10) fOverThreshold10 = new AliTPCCalPad("OverThreshold10","OverThreshold10");
580   if (!fOverThreshold20) fOverThreshold20 = new AliTPCCalPad("OverThreshold20","OverThreshold20");
581   if (!fOverThreshold30) fOverThreshold30 = new AliTPCCalPad("OverThreshold30","OverThreshold30");
582   if (!fHistQVsTimeSideA) {
583     fHistQVsTimeSideA  = new TProfile("hQVsTimeSideA", "Q vs time (side A); Time [Timebin]; Q [ADC ch]", 100, 0, 1000);
584     fHistQVsTimeSideA->SetDirectory(0);
585   }
586   if (!fHistQVsTimeSideC) {
587     fHistQVsTimeSideC  = new TProfile("hQVsTimeSideC", "Q vs time (side C); Time [Timebin]; Q [ADC ch]", 100, 0, 1000);
588     fHistQVsTimeSideC->SetDirectory(0);
589   }
590   if (!fHistQMaxVsTimeSideA) {
591     fHistQMaxVsTimeSideA  = new TProfile("hQMaxVsTimeSideA", "Q_{MAX} vs time (side A); Time [Timebin]; Q_{MAX} [ADC ch]", 100, 0, 1000);
592     fHistQMaxVsTimeSideA->SetDirectory(0);
593   }
594   if (!fHistQMaxVsTimeSideC) {
595     fHistQMaxVsTimeSideC  = new TProfile("hQMaxVsTimeSideC", "Q_{MAX} vs time (side C); Time [Timebin]; Q_{MAX} [ADC ch]", 100, 0, 1000);
596     fHistQMaxVsTimeSideC->SetDirectory(0);
597   }
598
599   // Make the arrays for expanding the data
600
601   if (!fAllBins)
602     MakeArrays();
603
604   //
605   // If Analyse has been previously called we need now to denormalize the data
606   // as more data is coming
607   //
608   if(fIsAnalysed == kTRUE) {
609     
610     const Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
611     const Float_t denormalization = Float_t(fEventCounter * nTimeBins);
612     fNoThreshold->Multiply(denormalization);  
613     
614     fMeanCharge->Multiply(fNLocalMaxima);
615     fMaxCharge->Multiply(fNLocalMaxima);
616     fNTimeBins->Multiply(fNLocalMaxima);
617     fNPads->Multiply(fNLocalMaxima);
618     fTimePosition->Multiply(fNLocalMaxima);
619     fIsAnalysed = kFALSE;
620   }
621
622   //
623   // TimeBin cut
624   //
625   if (iTimeBin<fFirstTimeBin) return 0;
626   if (iTimeBin>fLastTimeBin) return 0;
627   
628   // if pedestal calibrations are loaded subtract pedestals
629   if(fPedestal) {
630
631     Float_t ped = fPedestal->GetCalROC(iSector)->GetValue(iRow, iPad);
632     // Don't use data from pads where pedestals are abnormally small or large
633     if(ped<10 || ped>90)
634       return 0;
635     signal -= ped;
636   }
637   
638   // In fNoThreshold we fill all data to estimate the ZS volume
639   Float_t count = fNoThreshold->GetCalROC(iSector)->GetValue(iRow, iPad);
640   fNoThreshold->GetCalROC(iSector)->SetValue(iRow, iPad,count+1);
641   
642   // Require at least 3 ADC channels
643   if (signal < 3.0)
644     return 0;
645
646   // if noise calibrations are loaded require at least 3*sigmaNoise
647   if(fNoise) {
648     
649     Float_t noise = fNoise->GetCalROC(iSector)->GetValue(iRow, iPad);
650     
651     if(signal < noise*3.0)
652       return 0;
653   }
654
655   //
656   // This signal is ok and we store it in the cluster map
657   //
658
659   SetExpandDigit(iRow, iPad, iTimeBin, signal);
660   
661   return 1; // signal was accepted
662 }
663
664 //_____________________________________________________________________
665 void AliTPCdataQA::FindLocalMaxima(const Int_t iSector)
666 {
667   //
668   // This method is called after the data from each sector has been
669   // exapanded into an array
670   // Loop over the signals and identify local maxima and fill the
671   // calibration objects with the information
672   //
673
674   Int_t nLocalMaxima = 0;
675   const Int_t maxTimeBin = fTimeBinsMax+4; // Used to step between neighboring pads 
676                                            // Because we have tha pad-time data in a 
677                                            // 1d array
678
679   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
680
681     Float_t* allBins = fAllBins[iRow];
682     Int_t* sigBins   = fAllSigBins[iRow];
683     const Int_t nSigBins   = fAllNSigBins[iRow];
684     
685     for (Int_t iSig = 0; iSig < nSigBins; iSig++) {
686
687       Int_t bin  = sigBins[iSig];
688       Float_t *b = &allBins[bin];
689
690       //
691       // Now we check if this is a local maximum
692       //
693
694       Float_t qMax = b[0];
695
696       // First check that the charge is bigger than the threshold
697       if (qMax<5) 
698         continue;
699       
700       // Require at least one neighboring pad with signal
701       if (b[-maxTimeBin]+b[maxTimeBin]<=0) continue;
702
703       // Require at least one neighboring time bin with signal
704       if (b[-1]+b[1]<=0) continue;
705       
706       //
707       // Check that this is a local maximum
708       // Note that the checking is done so that if 2 charges has the same
709       // qMax then only 1 cluster is generated 
710       // (that is why there is BOTH > and >=)
711       //
712       if (b[-maxTimeBin]   >= qMax) continue;
713       if (b[-1  ]          >= qMax) continue; 
714       if (b[+maxTimeBin]   > qMax)  continue; 
715       if (b[+1  ]          > qMax)  continue; 
716       if (b[-maxTimeBin-1] >= qMax) continue;
717       if (b[+maxTimeBin-1] >= qMax) continue; 
718       if (b[+maxTimeBin+1] > qMax)  continue; 
719       if (b[-maxTimeBin+1] >= qMax) continue;
720       
721       //
722       // Now we accept the local maximum and fill the calibration/data objects
723       //
724       nLocalMaxima++;
725
726       Int_t iPad, iTimeBin;
727       GetPadAndTimeBin(bin, iPad, iTimeBin);
728       
729       Float_t count = fNLocalMaxima->GetCalROC(iSector)->GetValue(iRow, iPad);
730       fNLocalMaxima->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
731
732       count = fTimePosition->GetCalROC(iSector)->GetValue(iRow, iPad);
733       fTimePosition->GetCalROC(iSector)->SetValue(iRow, iPad, count+iTimeBin);
734       
735       Float_t charge = fMaxCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
736       fMaxCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qMax);
737       
738       if(qMax>=10) {
739         count = fOverThreshold10->GetCalROC(iSector)->GetValue(iRow, iPad);
740         fOverThreshold10->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
741       }
742       if(qMax>=20) {
743         count = fOverThreshold20->GetCalROC(iSector)->GetValue(iRow, iPad);
744         fOverThreshold20->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
745       }
746       if(qMax>=30) {
747         count = fOverThreshold30->GetCalROC(iSector)->GetValue(iRow, iPad);
748         fOverThreshold30->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
749       }
750
751       //
752       // Calculate the total charge as the sum over the region:
753       //
754       //    o o o o o
755       //    o i i i o
756       //    o i C i o
757       //    o i i i o
758       //    o o o o o
759       //
760       // with qmax at the center C.
761       //
762       // The inner charge (i) we always add, but we only add the outer
763       // charge (o) if the neighboring inner bin (i) has a signal.
764       //
765       Int_t minP = 0, maxP = 0, minT = 0, maxT = 0;
766       Float_t qTot = qMax;
767       for(Int_t i = -1; i<=1; i++) {
768         for(Int_t j = -1; j<=1; j++) {
769           
770           if(i==0 && j==0)
771             continue;
772           
773           Float_t charge1 = GetQ(b, i, j, maxTimeBin, minT, maxT, minP, maxP);
774           qTot += charge1;
775           if(charge1>0) {
776             // see if the next neighbor is also above threshold
777             if(i*j==0) {
778               qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP); 
779             } else {
780               // we are in a diagonal corner
781               qTot += GetQ(b,   i, 2*j, maxTimeBin, minT, maxT, minP, maxP); 
782               qTot += GetQ(b, 2*i,   j, maxTimeBin, minT, maxT, minP, maxP); 
783               qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP); 
784             }
785           }
786         }
787       }
788       
789       charge = fMeanCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
790       fMeanCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qTot);
791       
792       count = fNTimeBins->GetCalROC(iSector)->GetValue(iRow, iPad);
793       fNTimeBins->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxT-minT+1);
794       
795       count = fNPads->GetCalROC(iSector)->GetValue(iRow, iPad);
796       fNPads->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxP-minP+1);
797       
798       if((iSector%36)<18) { // A side
799         fHistQVsTimeSideA->Fill(iTimeBin, qTot);
800         fHistQMaxVsTimeSideA->Fill(iTimeBin, qMax);
801       } else {
802         fHistQVsTimeSideC->Fill(iTimeBin, qTot);
803         fHistQMaxVsTimeSideC->Fill(iTimeBin, qMax);      
804       }
805     } // end loop over signals
806   } // end loop over rows
807   
808   //  cout << "Number of local maximas found: " << nLocalMaxima << endl;
809 }
810
811 //_____________________________________________________________________
812 void AliTPCdataQA::Analyse()
813 {
814   //
815   //  Calculate calibration constants
816   //
817   
818   cout << "Analyse called" << endl;
819
820   if(fIsAnalysed == kTRUE) {
821     
822     cout << "No new data since Analyse was called last time" << endl;
823     return;
824   }
825
826   if(fEventCounter==0) {
827     
828     cout << "EventCounter == 0, Cannot analyse" << endl;
829     return;
830   }
831   
832   Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
833   cout << "EventCounter: " << fEventCounter << endl
834        << "TimeBins: " << nTimeBins << endl;
835
836   Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins);
837   fNoThreshold->Multiply(normalization);  
838   
839   fMeanCharge->Divide(fNLocalMaxima);
840   fMaxCharge->Divide(fNLocalMaxima);
841   fNTimeBins->Divide(fNLocalMaxima);
842   fNPads->Divide(fNLocalMaxima);
843   fTimePosition->Divide(fNLocalMaxima);
844
845   fIsAnalysed = kTRUE;
846 }
847
848
849 //_____________________________________________________________________
850 void AliTPCdataQA::MakeTree(const char *fname){
851   //
852   // Export result to the tree -located in the file
853   // This file can be analyzed using AliTPCCalibViewer
854   // 
855   AliTPCPreprocessorOnline preprocesor;
856
857   if (fNLocalMaxima) preprocesor.AddComponent(fNLocalMaxima);
858   if (fMaxCharge) preprocesor.AddComponent(fMaxCharge);  
859   if (fMeanCharge) preprocesor.AddComponent(fMeanCharge);  
860   if (fNoThreshold) preprocesor.AddComponent(fNoThreshold);
861   if (fNTimeBins) preprocesor.AddComponent(fNTimeBins);
862   if (fNPads) preprocesor.AddComponent(fNPads);
863   if (fTimePosition) preprocesor.AddComponent(fTimePosition);
864   if (fOverThreshold10) preprocesor.AddComponent(fOverThreshold10);
865   if (fOverThreshold20) preprocesor.AddComponent(fOverThreshold20);
866   if (fOverThreshold30) preprocesor.AddComponent(fOverThreshold30);
867
868   preprocesor.DumpToFile(fname);  
869 }
870
871
872 //_____________________________________________________________________
873 void AliTPCdataQA::MakeArrays(){
874   //
875   // The arrays for expanding the raw data are defined and 
876   // som parameters are intialised
877   //
878   AliTPCROC * roc = AliTPCROC::Instance();
879   //
880   // To make the array big enough for all sectors we take 
881   // the dimensions from the outer row of an OROC (the last sector)
882   //
883   fRowsMax     = roc->GetNRows(roc->GetNSector()-1);
884   fPadsMax     = roc->GetNPads(roc->GetNSector()-1,fRowsMax-1);
885   fTimeBinsMax = fLastTimeBin - fFirstTimeBin +1; 
886
887   //
888   // Since we have added 2 pads (TimeBins) before and after the real pads (TimeBins) 
889   // to make sure that we can always query the exanded table even when the 
890   // max is on the edge
891   //
892
893  
894   fAllBins = new Float_t*[fRowsMax];
895   fAllSigBins = new Int_t*[fRowsMax];
896   fAllNSigBins = new Int_t[fRowsMax];
897
898   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
899     //
900     Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4);  
901     fAllBins[iRow] = new Float_t[maxBin];
902     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); // set all values to 0
903     fAllSigBins[iRow] = new Int_t[maxBin];
904     fAllNSigBins[iRow] = 0;
905   }
906 }
907
908
909 //_____________________________________________________________________
910 void AliTPCdataQA::CleanArrays(){
911   //
912   //
913   //
914
915   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
916
917     // To speed up the performance by a factor 2 on cosmic data (and
918     // presumably pp data as well) where the ocupancy is low, the
919     // memset is only called if there is more than 1000 signals for a
920     // row (of the order 1% occupancy)
921     if(fAllNSigBins[iRow]<1000) {
922       
923       Float_t* allBins = fAllBins[iRow];
924       Int_t* sigBins   = fAllSigBins[iRow];
925       const Int_t nSignals = fAllNSigBins[iRow];
926       for(Int_t i = 0; i < nSignals; i++)
927         allBins[sigBins[i]]=0;      
928     } else {
929
930       Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4); 
931       memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
932     }
933
934     fAllNSigBins[iRow]=0;
935   }
936 }
937
938 //_____________________________________________________________________
939 void AliTPCdataQA::GetPadAndTimeBin(Int_t bin, Int_t& iPad, Int_t& iTimeBin){
940   //
941   // Return pad and timebin for a given bin
942   //
943   
944   //  Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
945   iTimeBin  = bin%(fTimeBinsMax+4);
946   iPad      = (bin-iTimeBin)/(fTimeBinsMax+4);
947
948   iPad     -= 2;
949   iTimeBin -= 2;
950   iTimeBin += fFirstTimeBin;
951
952   R__ASSERT(iPad>=0 && iPad<=fPadsMax);
953   R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
954 }
955
956 //_____________________________________________________________________
957 void AliTPCdataQA::SetExpandDigit(const Int_t iRow, Int_t iPad, 
958                                   Int_t iTimeBin, const Float_t signal) 
959 {
960   //
961   // 
962   //
963   R__ASSERT(iRow>=0 && iRow<fRowsMax);
964   R__ASSERT(iPad>=0 && iPad<=fPadsMax);
965   R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
966   
967   iTimeBin -= fFirstTimeBin;
968   iPad     += 2;
969   iTimeBin += 2;
970   
971   Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
972   fAllBins[iRow][bin] = signal;
973   fAllSigBins[iRow][fAllNSigBins[iRow]] = bin;
974   fAllNSigBins[iRow]++;
975 }
976
977 Float_t AliTPCdataQA::GetQ(const Float_t* adcArray, const Int_t time, 
978                            const Int_t pad, const Int_t maxTimeBins, 
979                            Int_t& timeMin, Int_t& timeMax, 
980                            Int_t& padMin,  Int_t& padMax) 
981 {
982   //
983   // This methods return the charge in the bin time+pad*maxTimeBins
984   // If the charge is above 0 it also updates the padMin, padMax, timeMin
985   // and timeMax if necessary
986   //
987   Float_t charge = adcArray[time + pad*maxTimeBins];
988   if(charge > 0) {
989     timeMin = TMath::Min(time, timeMin); timeMax = TMath::Max(time, timeMax);
990     padMin = TMath::Min(pad, padMin); padMax = TMath::Max(pad, padMax);
991   }
992   return charge; 
993 }
994
995 //______________________________________________________________________________
996 void AliTPCdataQA::Streamer(TBuffer &R__b)
997 {
998   // Automatic schema evolution was first used from revision 4
999   // Code based on:
1000   // http://root.cern.ch/root/roottalk/roottalk02/3207.html
1001
1002    UInt_t R__s, R__c;
1003    if (R__b.IsReading()) {
1004       Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1005       //we use the automatic algorithm for class version > 3
1006       if (R__v > 3) {
1007         AliTPCdataQA::Class()->ReadBuffer(R__b, this, R__v, R__s,
1008                                           R__c);
1009         return;
1010       }
1011       TH1F::Streamer(R__b);
1012       R__b >> fFirstTimeBin;
1013       R__b >> fLastTimeBin;
1014       R__b >> fAdcMin;
1015       R__b >> fAdcMax;
1016       R__b >> fNLocalMaxima;
1017       R__b >> fMaxCharge;
1018       R__b >> fMeanCharge;
1019       R__b >> fNoThreshold;
1020       R__b >> fNTimeBins;
1021       R__b >> fNPads;
1022       R__b >> fTimePosition;
1023       R__b >> fEventCounter;
1024       R__b >> fIsAnalysed;
1025       R__b.CheckByteCount(R__s, R__c, AliTPCdataQA::IsA());
1026    } else {
1027      AliTPCdataQA::Class()->WriteBuffer(R__b,this);
1028    }
1029 }