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