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