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