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