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
  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   fEventCounter(0),
 116   fIsAnalysed(kFALSE),
 117   fAllBins(0),
 118   fAllSigBins(0),
 119   fAllNSigBins(0),
 120   fRowsMax(0),
 121   fPadsMax(0),
 122   fTimeBinsMax(0)
 123 {
 124   //
 125   // default constructor
 126   //
 127 }
 128
 129 //_____________________________________________________________________
 130 AliTPCdataQA::AliTPCdataQA(AliRecoParam::EventSpecie_t es) :
 131 fFirstTimeBin(60),
 132 fLastTimeBin(1000),
 133 fAdcMin(1),
 134 fAdcMax(100),
 135 fMapping(NULL),
 136 fPedestal(0),
 137 fNoise(0),
 138 fNLocalMaxima(0),
 139 fMaxCharge(0),
 140 fMeanCharge(0),
 141 fNoThreshold(0),
 142 fNTimeBins(0),
 143 fNPads(0),
 144 fTimePosition(0),
 145 fOverThreshold10(0),
 146 fOverThreshold20(0),
 147 fOverThreshold30(0),
 148 fEventCounter(0),
 149 fIsAnalysed(kFALSE),
 150 fAllBins(0),
 151 fAllSigBins(0),
 152 fAllNSigBins(0),
 153 fRowsMax(0),
 154 fPadsMax(0),
 155 fTimeBinsMax(0)
 156 {
 157 // ctor creating the histogram
 158   char *  name = Form("TPCRAW_%s", AliRecoParam::GetEventSpecieName(es)) ;
 159   TH1F(name, name,100,0,100) ;
 160 }
 161
 162 //_____________________________________________________________________
 163 AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/
 164   TH1F(ped),
 165   fFirstTimeBin(ped.GetFirstTimeBin()),
 166   fLastTimeBin(ped.GetLastTimeBin()),
 167   fAdcMin(ped.GetAdcMin()),
 168   fAdcMax(ped.GetAdcMax()),
 169   fMapping(NULL),
 170   fPedestal(0),
 171   fNoise(0),
 172   fNLocalMaxima(0),
 173   fMaxCharge(0),
 174   fMeanCharge(0),
 175   fNoThreshold(0),
 176   fNTimeBins(0),
 177   fNPads(0),
 178   fTimePosition(0),
 179   fOverThreshold10(0),
 180   fOverThreshold20(0),
 181   fOverThreshold30(0),
 182   fEventCounter(ped.GetEventCounter()),
 183   fIsAnalysed(ped.GetIsAnalysed()),
 184   fAllBins(0),
 185   fAllSigBins(0),
 186   fAllNSigBins(0),
 187   fRowsMax(0),
 188   fPadsMax(0),
 189   fTimeBinsMax(0)
 190 {
 191   //
 192   // copy constructor
 193   //
 194   if(ped.GetNLocalMaxima())
 195     fNLocalMaxima  = new AliTPCCalPad(*ped.GetNLocalMaxima());
 196   if(ped.GetMaxCharge())
 197     fMaxCharge      = new AliTPCCalPad(*ped.GetMaxCharge());
 198   if(ped.GetMeanCharge())
 199     fMeanCharge     = new AliTPCCalPad(*ped.GetMeanCharge());
 200   if(ped.GetNoThreshold())
 201     fNoThreshold  = new AliTPCCalPad(*ped.GetNoThreshold());
 202   if(ped.GetNTimeBins())
 203     fNTimeBins  = new AliTPCCalPad(*ped.GetNTimeBins());
 204   if(ped.GetNPads())
 205     fNPads  = new AliTPCCalPad(*ped.GetNPads());
 206   if(ped.GetTimePosition())
 207     fTimePosition  = new AliTPCCalPad(*ped.GetTimePosition());
 208   if(ped.GetOverThreshold10())
 209     fOverThreshold10  = new AliTPCCalPad(*ped.GetOverThreshold10());
 210   if(ped.GetOverThreshold20())
 211     fOverThreshold20  = new AliTPCCalPad(*ped.GetOverThreshold20());
 212   if(ped.GetOverThreshold30())
 213     fOverThreshold30  = new AliTPCCalPad(*ped.GetOverThreshold30());
 214 }
 215
 216 //_____________________________________________________________________
 217 AliTPCdataQA::AliTPCdataQA(const TMap *config) : /*FOLD00*/
 218   TH1F("TPCRAW","TPCRAW",100,0,100),
 219   fFirstTimeBin(60),
 220   fLastTimeBin(1000),
 221   fAdcMin(1),
 222   fAdcMax(100),
 223   fMapping(NULL),
 224   fPedestal(0),
 225   fNoise(0),
 226   fNLocalMaxima(0),
 227   fMaxCharge(0),
 228   fMeanCharge(0),
 229   fNoThreshold(0),
 230   fNTimeBins(0),
 231   fNPads(0),
 232   fTimePosition(0),
 233   fOverThreshold10(0),
 234   fOverThreshold20(0),
 235   fOverThreshold30(0),
 236   fEventCounter(0),
 237   fIsAnalysed(kFALSE),
 238   fAllBins(0),
 239   fAllSigBins(0),
 240   fAllNSigBins(0),
 241   fRowsMax(0),
 242   fPadsMax(0),
 243   fTimeBinsMax(0)
 244 {
 245   //
 246   // default constructor
 247   //
 248   if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
 249   if (config->GetValue("LastTimeBin"))  fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
 250   if (config->GetValue("AdcMin"))       fAdcMin = ((TObjString*)config->GetValue("AdcMin"))->GetString().Atoi();
 251   if (config->GetValue("AdcMax"))       fAdcMax = ((TObjString*)config->GetValue("AdcMax"))->GetString().Atoi();
 252 }
 253
 254 //_____________________________________________________________________
 255 AliTPCdataQA& AliTPCdataQA::operator = (const  AliTPCdataQA &source)
 256 {
 257   //
 258   // assignment operator
 259   //
 260   if (&source == this) return *this;
 261   new (this) AliTPCdataQA(source);
 262
 263   return *this;
 264 }
 265
 266
 267 //_____________________________________________________________________
 268 AliTPCdataQA::~AliTPCdataQA() /*FOLD00*/
 269 {
 270   //
 271   // destructor
 272   //
 273
 274   // do not delete fMapping, because we do not own it.
 275   // do not delete fMapping, because we do not own it.
 276   // do not delete fNoise and fPedestal, because we do not own them.
 277
 278   delete fNLocalMaxima;
 279   delete fMaxCharge;
 280   delete fMeanCharge;
 281   delete fNoThreshold;
 282   delete fNTimeBins;
 283   delete fNPads;
 284   delete fTimePosition;
 285   delete fOverThreshold10;
 286   delete fOverThreshold20;
 287   delete fOverThreshold30;
 288
 289   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
 290     delete [] fAllBins[iRow];
 291     delete [] fAllSigBins[iRow];
 292   }
 293   delete [] fAllBins;
 294   delete [] fAllSigBins;
 295   delete [] fAllNSigBins;
 296 }
 297 //_____________________________________________________________________
 298 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStreamV3 *rawStreamV3)
 299 {
 300   //
 301   // Event Processing loop - AliTPCRawStreamV3
 302   //
 303   Bool_t withInput = kFALSE;
 304   Int_t nSignals = 0;
 305   Int_t lastSector = -1;
 306
 307   while ( rawStreamV3->NextDDL() ){
 308     while ( rawStreamV3->NextChannel() ){
 309       Int_t iSector = rawStreamV3->GetSector(); //  current sector
 310       Int_t iRow    = rawStreamV3->GetRow();    //  current row
 311       Int_t iPad    = rawStreamV3->GetPad();    //  current pad
 312       if (iRow<0 || iPad<0) continue;
 313       // Call local maxima finder if the data is in a new sector
 314       if(iSector != lastSector) {
 315
 316         if(nSignals>0)
 317           FindLocalMaxima(lastSector);
 318
 319         CleanArrays();
 320         lastSector = iSector;
 321         nSignals = 0;
 322       }
 323
 324       while ( rawStreamV3->NextBunch() ){
 325         Int_t  startTbin    = (Int_t)rawStreamV3->GetStartTimeBin();
 326         Int_t  bunchlength  = (Int_t)rawStreamV3->GetBunchLength();
 327         const UShort_t *sig = rawStreamV3->GetSignals();
 328
 329         for (Int_t iTimeBin = 0; iTimeBin<bunchlength; iTimeBin++){
 330           Float_t signal=(Float_t)sig[iTimeBin];
 331           nSignals += Update(iSector,iRow,iPad,startTbin--,signal);
 332           withInput = kTRUE;
 333         }
 334       }
 335     }
 336   }
 337
 338   if (lastSector>=0&&nSignals>0)
 339     FindLocalMaxima(lastSector);
 340
 341   return withInput;
 342 }
 343
 344 //_____________________________________________________________________
 345 Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader)
 346 {
 347   //
 348   //  Event processing loop - AliRawReader
 349   //
 350   AliTPCRawStreamV3 *rawStreamV3 = new AliTPCRawStreamV3(rawReader, (AliAltroMapping**)fMapping);
 351   Bool_t res=ProcessEvent(rawStreamV3);
 352   delete rawStreamV3;
 353   if(res)
 354     fEventCounter++; // only increment event counter if there is TPC data
 355   return res;
 356 }
 357
 358 //_____________________________________________________________________
 359 Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
 360 {
 361   //
 362   // Event Processing loop - AliTPCRawStream
 363   //
 364   Bool_t withInput = kFALSE;
 365   Int_t nSignals = 0;
 366   Int_t lastSector = -1;
 367
 368   while ( rawStreamFast->NextDDL() ){
 369     while ( rawStreamFast->NextChannel() ){
 370
 371       Int_t iSector  = rawStreamFast->GetSector(); //  current sector
 372       Int_t iRow     = rawStreamFast->GetRow();    //  current row
 373       Int_t iPad     = rawStreamFast->GetPad();    //  current pad
 374   // Call local maxima finder if the data is in a new sector
 375       if(iSector != lastSector) {
 376
 377         if(nSignals>0)
 378           FindLocalMaxima(lastSector);
 379
 380         CleanArrays();
 381         lastSector = iSector;
 382         nSignals = 0;
 383       }
 384
 385       while ( rawStreamFast->NextBunch() ){
 386         Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
 387         Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
 388
 389         for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
 390           Float_t signal = rawStreamFast->GetSignals()[iTimeBin-startTbin];
 391           nSignals += Update(iSector,iRow,iPad,iTimeBin+1,signal);
 392           withInput = kTRUE;
 393         }
 394       }
 395     }
 396   }
 397
 398   return withInput;
 399 }
 400 //_____________________________________________________________________
 401 Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *rawReader)
 402 {
 403   //
 404   //  Event processing loop - AliRawReader
 405   //
 406   AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
 407   Bool_t res=ProcessEventFast(rawStreamFast);
 408   if(res)
 409     fEventCounter++; // only increment event counter if there is TPC data
 410                      // otherwise Analyse (called in QA) fails
 411
 412   delete rawStreamFast;
 413   return res;
 414 }
 415
 416 //_____________________________________________________________________
 417 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream)
 418 {
 419   //
 420   // Event Processing loop - AliTPCRawStream
 421   //
 422
 423   Bool_t withInput = kFALSE;
 424   Int_t nSignals = 0;
 425   Int_t lastSector = -1;
 426
 427   while (rawStream->Next()) {
 428
 429     Int_t iSector  = rawStream->GetSector();      //  current ROC
 430     Int_t iRow     = rawStream->GetRow();         //  current row
 431     Int_t iPad     = rawStream->GetPad();         //  current pad
 432     Int_t iTimeBin = rawStream->GetTime();        //  current time bin
 433     Float_t signal = rawStream->GetSignal();      //  current ADC signal
 434
 435     // Call local maxima finder if the data is in a new sector
 436     if(iSector != lastSector) {
 437
 438       if(nSignals>0)
 439         FindLocalMaxima(lastSector);
 440
 441       CleanArrays();
 442       lastSector = iSector;
 443       nSignals = 0;
 444     }
 445
 446     // Sometimes iRow==-1 if there is problems to read the data
 447     if(iRow>=0) {
 448       nSignals += Update(iSector,iRow,iPad,iTimeBin,signal);
 449       withInput = kTRUE;
 450     }
 451   }
 452
 453   if (lastSector>=0&&nSignals>0)
 454     FindLocalMaxima(lastSector);
 455
 456   return withInput;
 457 }
 458
 459
 460 //_____________________________________________________________________
 461 Bool_t AliTPCdataQA::ProcessEventOld(AliRawReader *rawReader)
 462 {
 463   //
 464   //  Event processing loop - AliRawReader
 465   //
 466
 467   // if fMapping is NULL the rawstream will crate its own mapping
 468   AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
 469   rawReader->Select("TPC");
 470   Bool_t res =  ProcessEvent(&rawStream);
 471
 472   if(res)
 473     fEventCounter++; // only increment event counter if there is TPC data
 474                      // otherwise Analyse (called in QA) fails
 475   return res;
 476 }
 477
 478
 479 //_____________________________________________________________________
 480 Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *event)
 481 {
 482   //
 483   //  process date event
 484   //
 485
 486   AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 487   Bool_t result=ProcessEvent(rawReader);
 488   delete rawReader;
 489   return result;
 490 }
 491
 492
 493
 494 //_____________________________________________________________________
 495 void AliTPCdataQA::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
 496 {
 497   //
 498   //  Write class to file
 499   //
 500
 501   TString sDir(dir);
 502   TString option;
 503
 504   if ( append )
 505     option = "update";
 506   else
 507     option = "recreate";
 508
 509   TDirectory *backup = gDirectory;
 510   TFile f(filename,option.Data());
 511   f.cd();
 512   if ( !sDir.IsNull() ){
 513     f.mkdir(sDir.Data());
 514     f.cd(sDir);
 515   }
 516   this->Write();
 517   f.Close();
 518
 519   if ( backup ) backup->cd();
 520 }
 521
 522
 523 //_____________________________________________________________________
 524 Int_t AliTPCdataQA::Update(const Int_t iSector, /*FOLD00*/
 525                            const Int_t iRow,
 526                            const Int_t iPad,
 527                            const Int_t iTimeBin,
 528                            Float_t signal)
 529 {
 530   //
 531   // Signal filling method
 532   //
 533
 534   //
 535   // Define the calibration objects the first time Update is called
 536   // NB! This has to be done first even if the data is rejected by the time
 537   // cut to make sure that the objects are available in Analyse
 538   //
 539   if (!fNLocalMaxima) fNLocalMaxima = new AliTPCCalPad("NLocalMaxima","NLocalMaxima");
 540   if (!fMaxCharge) fMaxCharge = new AliTPCCalPad("MaxCharge","MaxCharge");
 541   if (!fMeanCharge) fMeanCharge = new AliTPCCalPad("MeanCharge","MeanCharge");
 542   if (!fNoThreshold) fNoThreshold = new AliTPCCalPad("NoThreshold","NoThreshold");
 543   if (!fNTimeBins) fNTimeBins = new AliTPCCalPad("NTimeBins","NTimeBins");
 544   if (!fNPads) fNPads = new AliTPCCalPad("NPads","NPads");
 545   if (!fTimePosition) fTimePosition = new AliTPCCalPad("TimePosition","TimePosition");
 546   if (!fOverThreshold10) fOverThreshold10 = new AliTPCCalPad("OverThreshold10","OverThreshold10");
 547   if (!fOverThreshold20) fOverThreshold20 = new AliTPCCalPad("OverThreshold20","OverThreshold20");
 548   if (!fOverThreshold30) fOverThreshold30 = new AliTPCCalPad("OverThreshold30","OverThreshold30");
 549   // Make the arrays for expanding the data
 550
 551   if (!fAllBins)
 552     MakeArrays();
 553
 554   //
 555   // If Analyse has been previously called we need now to denormalize the data
 556   // as more data is coming
 557   //
 558   if(fIsAnalysed == kTRUE) {
 559
 560     const Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
 561     const Float_t denormalization = Float_t(fEventCounter * nTimeBins);
 562     fNoThreshold->Multiply(denormalization);
 563
 564     fMeanCharge->Multiply(fNLocalMaxima);
 565     fMaxCharge->Multiply(fNLocalMaxima);
 566     fNTimeBins->Multiply(fNLocalMaxima);
 567     fNPads->Multiply(fNLocalMaxima);
 568     fTimePosition->Multiply(fNLocalMaxima);
 569     fIsAnalysed = kFALSE;
 570   }
 571
 572   //
 573   // TimeBin cut
 574   //
 575   if (iTimeBin<fFirstTimeBin) return 0;
 576   if (iTimeBin>fLastTimeBin) return 0;
 577
 578   // if pedestal calibrations are loaded subtract pedestals
 579   if(fPedestal) {
 580
 581     Float_t ped = fPedestal->GetCalROC(iSector)->GetValue(iRow, iPad);
 582     // Don't use data from pads where pedestals are abnormally small or large
 583     if(ped<10 || ped>90)
 584       return 0;
 585     signal -= ped;
 586   }
 587
 588   // In fNoThreshold we fill all data to estimate the ZS volume
 589   Float_t count = fNoThreshold->GetCalROC(iSector)->GetValue(iRow, iPad);
 590   fNoThreshold->GetCalROC(iSector)->SetValue(iRow, iPad,count+1);
 591
 592   // Require at least 3 ADC channels
 593   if (signal < 3.0)
 594     return 0;
 595
 596   // if noise calibrations are loaded require at least 3*sigmaNoise
 597   if(fNoise) {
 598
 599     Float_t noise = fNoise->GetCalROC(iSector)->GetValue(iRow, iPad);
 600
 601     if(signal < noise*3.0)
 602       return 0;
 603   }
 604
 605   //
 606   // This signal is ok and we store it in the cluster map
 607   //
 608
 609   SetExpandDigit(iRow, iPad, iTimeBin, signal);
 610
 611   return 1; // signal was accepted
 612 }
 613
 614 //_____________________________________________________________________
 615 void AliTPCdataQA::FindLocalMaxima(const Int_t iSector)
 616 {
 617   //
 618   // This method is called after the data from each sector has been
 619   // exapanded into an array
 620   // Loop over the signals and identify local maxima and fill the
 621   // calibration objects with the information
 622   //
 623
 624   Int_t nLocalMaxima = 0;
 625   const Int_t maxTimeBin = fTimeBinsMax+4; // Used to step between neighboring pads
 626                                            // Because we have tha pad-time data in a
 627                                            // 1d array
 628
 629   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
 630
 631     Float_t* allBins = fAllBins[iRow];
 632     Int_t* sigBins   = fAllSigBins[iRow];
 633     const Int_t nSigBins   = fAllNSigBins[iRow];
 634
 635     for (Int_t iSig = 0; iSig < nSigBins; iSig++) {
 636
 637       Int_t bin  = sigBins[iSig];
 638       Float_t *b = &allBins[bin];
 639
 640       //
 641       // Now we check if this is a local maximum
 642       //
 643
 644       Float_t qMax = b[0];
 645
 646       // First check that the charge is bigger than the threshold
 647       if (qMax<5)
 648         continue;
 649
 650       // Require at least one neighboring pad with signal
 651       if (b[-maxTimeBin]+b[maxTimeBin]<=0) continue;
 652
 653       // Require at least one neighboring time bin with signal
 654       if (b[-1]+b[1]<=0) continue;
 655
 656       //
 657       // Check that this is a local maximum
 658       // Note that the checking is done so that if 2 charges has the same
 659       // qMax then only 1 cluster is generated
 660       // (that is why there is BOTH > and >=)
 661       //
 662       if (b[-maxTimeBin]   >= qMax) continue;
 663       if (b[-1  ]          >= qMax) continue;
 664       if (b[+maxTimeBin]   > qMax)  continue;
 665       if (b[+1  ]          > qMax)  continue;
 666       if (b[-maxTimeBin-1] >= qMax) continue;
 667       if (b[+maxTimeBin-1] >= qMax) continue;
 668       if (b[+maxTimeBin+1] > qMax)  continue;
 669       if (b[-maxTimeBin+1] >= qMax) continue;
 670
 671       //
 672       // Now we accept the local maximum and fill the calibration/data objects
 673       //
 674       nLocalMaxima++;
 675
 676       Int_t iPad, iTimeBin;
 677       GetPadAndTimeBin(bin, iPad, iTimeBin);
 678
 679       Float_t count = fNLocalMaxima->GetCalROC(iSector)->GetValue(iRow, iPad);
 680       fNLocalMaxima->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
 681
 682       count = fTimePosition->GetCalROC(iSector)->GetValue(iRow, iPad);
 683       fTimePosition->GetCalROC(iSector)->SetValue(iRow, iPad, count+iTimeBin);
 684
 685       Float_t charge = fMaxCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
 686       fMaxCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qMax);
 687
 688       if(qMax>=10) {
 689         count = fOverThreshold10->GetCalROC(iSector)->GetValue(iRow, iPad);
 690         fOverThreshold10->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
 691       }
 692       if(qMax>=20) {
 693         count = fOverThreshold20->GetCalROC(iSector)->GetValue(iRow, iPad);
 694         fOverThreshold20->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
 695       }
 696       if(qMax>=30) {
 697         count = fOverThreshold30->GetCalROC(iSector)->GetValue(iRow, iPad);
 698         fOverThreshold30->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
 699       }
 700
 701       //
 702       // Calculate the total charge as the sum over the region:
 703       //
 704       //    o o o o o
 705       //    o i i i o
 706       //    o i C i o
 707       //    o i i i o
 708       //    o o o o o
 709       //
 710       // with qmax at the center C.
 711       //
 712       // The inner charge (i) we always add, but we only add the outer
 713       // charge (o) if the neighboring inner bin (i) has a signal.
 714       //
 715       Int_t minP = 0, maxP = 0, minT = 0, maxT = 0;
 716       Float_t qTot = qMax;
 717       for(Int_t i = -1; i<=1; i++) {
 718         for(Int_t j = -1; j<=1; j++) {
 719
 720           if(i==0 && j==0)
 721             continue;
 722
 723           Float_t charge1 = GetQ(b, i, j, maxTimeBin, minT, maxT, minP, maxP);
 724           qTot += charge1;
 725           if(charge1>0) {
 726             // see if the next neighbor is also above threshold
 727             if(i*j==0) {
 728               qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
 729             } else {
 730               // we are in a diagonal corner
 731               qTot += GetQ(b,   i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
 732               qTot += GetQ(b, 2*i,   j, maxTimeBin, minT, maxT, minP, maxP);
 733               qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
 734             }
 735           }
 736         }
 737       }
 738
 739       charge = fMeanCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
 740       fMeanCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qTot);
 741
 742       count = fNTimeBins->GetCalROC(iSector)->GetValue(iRow, iPad);
 743       fNTimeBins->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxT-minT+1);
 744
 745       count = fNPads->GetCalROC(iSector)->GetValue(iRow, iPad);
 746       fNPads->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxP-minP+1);
 747
 748     } // end loop over signals
 749   } // end loop over rows
 750
 751   //  cout << "Number of local maximas found: " << nLocalMaxima << endl;
 752 }
 753
 754 //_____________________________________________________________________
 755 void AliTPCdataQA::Analyse()
 756 {
 757   //
 758   //  Calculate calibration constants
 759   //
 760
 761   cout << "Analyse called" << endl;
 762
 763   if(fIsAnalysed == kTRUE) {
 764
 765     cout << "No new data since Analyse was called last time" << endl;
 766     return;
 767   }
 768
 769   if(fEventCounter==0) {
 770
 771     cout << "EventCounter == 0, Cannot analyse" << endl;
 772     return;
 773   }
 774
 775   Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
 776   cout << "EventCounter: " << fEventCounter << endl
 777        << "TimeBins: " << nTimeBins << endl;
 778
 779   Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins);
 780   fNoThreshold->Multiply(normalization);
 781
 782   fMeanCharge->Divide(fNLocalMaxima);
 783   fMaxCharge->Divide(fNLocalMaxima);
 784   fNTimeBins->Divide(fNLocalMaxima);
 785   fNPads->Divide(fNLocalMaxima);
 786   fTimePosition->Divide(fNLocalMaxima);
 787
 788   fIsAnalysed = kTRUE;
 789 }
 790
 791
 792 //_____________________________________________________________________
 793 void AliTPCdataQA::MakeTree(const char *fname){
 794   //
 795   // Export result to the tree -located in the file
 796   // This file can be analyzed using AliTPCCalibViewer
 797   //
 798   AliTPCPreprocessorOnline preprocesor;
 799
 800   if (fNLocalMaxima) preprocesor.AddComponent(fNLocalMaxima);
 801   if (fMaxCharge) preprocesor.AddComponent(fMaxCharge);
 802   if (fMeanCharge) preprocesor.AddComponent(fMeanCharge);
 803   if (fNoThreshold) preprocesor.AddComponent(fNoThreshold);
 804   if (fNTimeBins) preprocesor.AddComponent(fNTimeBins);
 805   if (fNPads) preprocesor.AddComponent(fNPads);
 806   if (fTimePosition) preprocesor.AddComponent(fTimePosition);
 807   if (fOverThreshold10) preprocesor.AddComponent(fOverThreshold10);
 808   if (fOverThreshold20) preprocesor.AddComponent(fOverThreshold20);
 809   if (fOverThreshold30) preprocesor.AddComponent(fOverThreshold30);
 810
 811   preprocesor.DumpToFile(fname);
 812 }
 813
 814
 815 //_____________________________________________________________________
 816 void AliTPCdataQA::MakeArrays(){
 817   //
 818   // The arrays for expanding the raw data are defined and
 819   // som parameters are intialised
 820   //
 821   AliTPCROC * roc = AliTPCROC::Instance();
 822   //
 823   // To make the array big enough for all sectors we take
 824   // the dimensions from the outer row of an OROC (the last sector)
 825   //
 826   fRowsMax     = roc->GetNRows(roc->GetNSector()-1);
 827   fPadsMax     = roc->GetNPads(roc->GetNSector()-1,fRowsMax-1);
 828   fTimeBinsMax = fLastTimeBin - fFirstTimeBin +1;
 829
 830   //
 831   // Since we have added 2 pads (TimeBins) before and after the real pads (TimeBins)
 832   // to make sure that we can always query the exanded table even when the
 833   // max is on the edge
 834   //
 835
 836
 837   fAllBins = new Float_t*[fRowsMax];
 838   fAllSigBins = new Int_t*[fRowsMax];
 839   fAllNSigBins = new Int_t[fRowsMax];
 840
 841   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
 842     //
 843     Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4);
 844     fAllBins[iRow] = new Float_t[maxBin];
 845     memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); // set all values to 0
 846     fAllSigBins[iRow] = new Int_t[maxBin];
 847     fAllNSigBins[iRow] = 0;
 848   }
 849 }
 850
 851
 852 //_____________________________________________________________________
 853 void AliTPCdataQA::CleanArrays(){
 854   //
 855   //
 856   //
 857
 858   for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
 859
 860     // To speed up the performance by a factor 2 on cosmic data (and
 861     // presumably pp data as well) where the ocupancy is low, the
 862     // memset is only called if there is more than 1000 signals for a
 863     // row (of the order 1% occupancy)
 864     if(fAllNSigBins[iRow]<1000) {
 865
 866       Float_t* allBins = fAllBins[iRow];
 867       Int_t* sigBins   = fAllSigBins[iRow];
 868       const Int_t nSignals = fAllNSigBins[iRow];
 869       for(Int_t i = 0; i < nSignals; i++)
 870         allBins[sigBins[i]]=0;
 871
 872     } else {
 873       Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4);
 874       memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
 875     }
 876
 877     fAllNSigBins[iRow]=0;
 878   }
 879 }
 880
 881 //_____________________________________________________________________
 882 void AliTPCdataQA::GetPadAndTimeBin(Int_t bin, Int_t& iPad, Int_t& iTimeBin){
 883   //
 884   // Return pad and timebin for a given bin
 885   //
 886
 887   //  Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
 888   iTimeBin  = bin%(fTimeBinsMax+4);
 889   iPad      = (bin-iTimeBin)/(fTimeBinsMax+4);
 890
 891   iPad     -= 2;
 892   iTimeBin -= 2;
 893   iTimeBin += fFirstTimeBin;
 894
 895   R__ASSERT(iPad>=0 && iPad<=fPadsMax);
 896   R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
 897 }
 898
 899 //_____________________________________________________________________
 900 void AliTPCdataQA::SetExpandDigit(const Int_t iRow, Int_t iPad,
 901                                   Int_t iTimeBin, const Float_t signal)
 902 {
 903   //
 904   //
 905   //
 906   R__ASSERT(iRow>=0 && iRow<fRowsMax);
 907   R__ASSERT(iPad>=0 && iPad<=fPadsMax);
 908   R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
 909
 910   iTimeBin -= fFirstTimeBin;
 911   iPad     += 2;
 912   iTimeBin += 2;
 913
 914   Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
 915   fAllBins[iRow][bin] = signal;
 916   fAllSigBins[iRow][fAllNSigBins[iRow]] = bin;
 917   fAllNSigBins[iRow]++;
 918 }
 919
 920 Float_t AliTPCdataQA::GetQ(const Float_t* adcArray, const Int_t time,
 921                            const Int_t pad, const Int_t maxTimeBins,
 922                            Int_t& timeMin, Int_t& timeMax,
 923                            Int_t& padMin,  Int_t& padMax)
 924 {
 925   //
 926   // This methods return the charge in the bin time+pad*maxTimeBins
 927   // If the charge is above 0 it also updates the padMin, padMax, timeMin
 928   // and timeMax if necessary
 929   //
 930   Float_t charge = adcArray[time + pad*maxTimeBins];
 931   if(charge > 0) {
 932     timeMin = TMath::Min(time, timeMin); timeMax = TMath::Max(time, timeMax);
 933     padMin = TMath::Min(pad, padMin); padMax = TMath::Max(pad, padMax);
 934   }
 935   return charge;
 936 }