TPC/AliTPCdataQA.cxx

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