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