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 }