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