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