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