TPC/AliTPCCalibPulser.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 /* $Id$ */
  17
  18 //Root includes
  19 #include <TH1F.h>
  20 #include <TH2S.h>
  21 #include <TString.h>
  22 #include <TVectorF.h>
  23 #include <TMath.h>
  24
  25 #include <TDirectory.h>
  26 #include <TSystem.h>
  27 #include <TFile.h>
  28
  29 //AliRoot includes
  30 #include "AliRawReader.h"
  31 #include "AliRawReaderRoot.h"
  32 #include "AliRawReaderDate.h"
  33 #include "AliTPCRawStream.h"
  34 #include "AliTPCCalROC.h"
  35 #include "AliTPCCalPad.h"
  36 #include "AliTPCROC.h"
  37 #include "AliTPCParam.h"
  38 #include "AliTPCCalibPulser.h"
  39 #include "AliTPCcalibDB.h"
  40 #include "AliMathBase.h"
  41 #include "TTreeStream.h"
  42
  43 //date
  44 #include "event.h"
  45
  46
  47 ///////////////////////////////////////////////////////////////////////////////////////
  48 //          Implementation of the TPC pulser calibration
  49 //
  50 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
  51 //
  52 //
  53 /***************************************************************************
  54  *                      Class Description                                  *
  55  ***************************************************************************
  56
  57  The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
  58  runs performed with the calibration pulser.
  59
  60  The information retrieved is
  61  - Time0 differences
  62  - Signal width differences
  63  - Amplification variations
  64
  65  the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
  66  one chip and somewhat large between different chips.
  67
  68  Histograms:
  69    For each ROC three TH2S histos 'Reference Histograms'  (ROC channel vs. [Time0, signal width, Q sum]) is created when
  70    it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
  71    TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
  72
  73
  74  Working principle:
  75  ------------------
  76  Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
  77  (see below). These in the end call the Update(...) function.
  78
  79  - the Update(...) function:
  80    In this function the array fPadSignal is filled with the adc signals between the specified range
  81    fFirstTimeBin and fLastTimeBin for the current pad.
  82    before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
  83    stored in fPadSignal.
  84
  85    - the ProcessPad() function:
  86      Find Pedestal and Noise information
  87      - use database information which has to be set by calling
  88        SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
  89      - if no information from the pedestal data base
  90        is available the informaion is calculated on the fly ( see FindPedestal() function )
  91
  92      Find the Pulser signal information
  93      - calculate  mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
  94        the Q sum is scaled by pad area
  95        (see FindPulserSignal(...) function)
  96
  97      Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
  98      Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
  99
 100  At the end of each event the EndEvent() function is called
 101
 102  - the EndEvent() function:
 103    calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
 104    This is done to overcome syncronisation problems between the trigger and the fec clock.
 105
 106  After accumulating the desired statistics the Analyse() function has to be called.
 107  - the Analyse() function
 108    Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
 109    the AliMathBase::GetCOG(...) function, and the calibration
 110    storage classes (AliTPCCalROC) are filled for each ROC.
 111    The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
 112    fCalRocArrayQ;
 113
 114
 115
 116  User interface for filling data:
 117  --------------------------------
 118
 119  To Fill information one of the following functions can be used:
 120
 121  Bool_t ProcessEvent(eventHeaderStruct *event);
 122    - process Date event
 123    - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
 124
 125  Bool_t ProcessEvent(AliRawReader *rawReader);
 126    - process AliRawReader event
 127    - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
 128
 129  Bool_t ProcessEvent(AliTPCRawStream *rawStream);
 130    - process event from AliTPCRawStream
 131    - call Update function for signal filling
 132
 133  Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
 134               iPad,  const Int_t iTimeBin, const Float_t signal);
 135    - directly  fill signal information (sector, row, pad, time bin, pad)
 136      to the reference histograms
 137
 138  It is also possible to merge two independently taken calibrations using the function
 139
 140  void Merge(AliTPCCalibPulser *sig)
 141    - copy histograms in 'sig' if the do not exist in this instance
 142    - Add histograms in 'sig' to the histograms in this instance if the allready exist
 143    - After merging call Analyse again!
 144
 145
 146
 147  -- example: filling data using root raw data:
 148  void fillSignal(Char_t *filename)
 149  {
 150     rawReader = new AliRawReaderRoot(fileName);
 151     if ( !rawReader ) return;
 152     AliTPCCalibPulser *calib = new AliTPCCalibPulser;
 153     while (rawReader->NextEvent()){
 154       calib->ProcessEvent(rawReader);
 155     }
 156     calib->Analyse();
 157     calib->DumpToFile("SignalData.root");
 158     delete rawReader;
 159     delete calib;
 160  }
 161
 162
 163  What kind of information is stored and how to retrieve them:
 164  ------------------------------------------------------------
 165
 166  - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
 167
 168    TH2F *GetHistoT0(Int_t sector);
 169    TH2F *GetHistoRMS(Int_t sector);
 170    TH2F *GetHistoQ(Int_t sector);
 171
 172  - Accessing the calibration storage objects:
 173
 174    AliTPCCalROC *GetCalRocT0(Int_t sector);   // for the Time0 values
 175    AliTPCCalROC *GetCalRocRMS(Int_t sector);  // for the signal width values
 176    AliTPCCalROC *GetCalRocQ(Int_t sector);    // for the Q sum values
 177
 178    example for visualisation:
 179    if the file "SignalData.root" was created using the above example one could do the following:
 180
 181    TFile fileSignal("SignalData.root")
 182    AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
 183    sig->GetCalRocT0(0)->Draw("colz");
 184    sig->GetCalRocRMS(0)->Draw("colz");
 185
 186    or use the AliTPCCalPad functionality:
 187    AliTPCCalPad padT0(ped->GetCalPadT0());
 188    AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
 189    padT0->MakeHisto2D()->Draw("colz");       //Draw A-Side Time0 Information
 190    padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
 191 */
 192
 193
 194 ClassImp(AliTPCCalibPulser) /*FOLD00*/
 195
 196 AliTPCCalibPulser::AliTPCCalibPulser() : /*FOLD00*/
 197     TObject(),
 198     fFirstTimeBin(60),
 199     fLastTimeBin(120),
 200     fNbinsT0(200),
 201     fXminT0(-2),
 202     fXmaxT0(2),
 203     fNbinsQ(200),
 204     fXminQ(1),
 205     fXmaxQ(40),
 206     fNbinsRMS(100),
 207     fXminRMS(0.1),
 208     fXmaxRMS(5.1),
 209     fLastSector(-1),
 210     fOldRCUformat(kTRUE),
 211     fROC(AliTPCROC::Instance()),
 212     fParam(new AliTPCParam),
 213     fPedestalTPC(0x0),
 214     fPadNoiseTPC(0x0),
 215     fPedestalROC(0x0),
 216     fPadNoiseROC(0x0),
 217     fCalRocArrayT0(72),
 218     fCalRocArrayQ(72),
 219     fCalRocArrayRMS(72),
 220     fCalRocArrayOutliers(72),
 221     fHistoQArray(72),
 222     fHistoT0Array(72),
 223     fHistoRMSArray(72),
 224     fPadTimesArrayEvent(72),
 225     fPadQArrayEvent(72),
 226     fPadRMSArrayEvent(72),
 227     fPadPedestalArrayEvent(72),
 228     fCurrentChannel(-1),
 229     fCurrentSector(-1),
 230     fCurrentRow(-1),
 231     fMaxPadSignal(-1),
 232     fMaxTimeBin(-1),
 233     fPadSignal(1024),
 234     fPadPedestal(0),
 235     fPadNoise(0),
 236     fVTime0Offset(72),
 237     fVTime0OffsetCounter(72),
 238     fEvent(-1),
 239     fDebugStreamer(0x0),
 240     fDebugLevel(0)
 241 {
 242     //
 243     // AliTPCSignal default constructor
 244     //
 245
 246 }
 247 //_____________________________________________________________________
 248 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
 249     TObject(sig),
 250     fFirstTimeBin(sig.fFirstTimeBin),
 251     fLastTimeBin(sig.fLastTimeBin),
 252     fNbinsT0(sig.fNbinsT0),
 253     fXminT0(sig.fXminT0),
 254     fXmaxT0(sig.fXmaxT0),
 255     fNbinsQ(sig.fNbinsQ),
 256     fXminQ(sig.fXminQ),
 257     fXmaxQ(sig.fXmaxQ),
 258     fNbinsRMS(sig.fNbinsRMS),
 259     fXminRMS(sig.fXminRMS),
 260     fXmaxRMS(sig.fXmaxRMS),
 261     fLastSector(-1),
 262     fOldRCUformat(kTRUE),
 263     fROC(AliTPCROC::Instance()),
 264     fParam(new AliTPCParam),
 265     fPedestalTPC(0x0),
 266     fPadNoiseTPC(0x0),
 267     fPedestalROC(0x0),
 268     fPadNoiseROC(0x0),
 269     fCalRocArrayT0(72),
 270     fCalRocArrayQ(72),
 271     fCalRocArrayRMS(72),
 272     fCalRocArrayOutliers(72),
 273     fHistoQArray(72),
 274     fHistoT0Array(72),
 275     fHistoRMSArray(72),
 276     fPadTimesArrayEvent(72),
 277     fPadQArrayEvent(72),
 278     fPadRMSArrayEvent(72),
 279     fPadPedestalArrayEvent(72),
 280     fCurrentChannel(-1),
 281     fCurrentSector(-1),
 282     fCurrentRow(-1),
 283     fMaxPadSignal(-1),
 284     fMaxTimeBin(-1),
 285     fPadSignal(1024),
 286     fPadPedestal(0),
 287     fPadNoise(0),
 288     fVTime0Offset(72),
 289     fVTime0OffsetCounter(72),
 290     fEvent(-1),
 291     fDebugStreamer(0x0),
 292     fDebugLevel(sig.fDebugLevel)
 293 {
 294     //
 295     // AliTPCSignal default constructor
 296     //
 297
 298     for (Int_t iSec = 0; iSec < 72; ++iSec){
 299         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 300         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 301         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 302         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 303
 304         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 305         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 306         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 307
 308         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 309         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 310         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 311         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 312
 313         if ( hQ != 0x0 ){
 314             TH2S *hNew = new TH2S(*hQ);
 315             hNew->SetDirectory(0);
 316             fHistoQArray.AddAt(hNew,iSec);
 317         }
 318         if ( hT0 != 0x0 ){
 319             TH2S *hNew = new TH2S(*hT0);
 320             hNew->SetDirectory(0);
 321             fHistoQArray.AddAt(hNew,iSec);
 322         }
 323         if ( hRMS != 0x0 ){
 324             TH2S *hNew = new TH2S(*hRMS);
 325             hNew->SetDirectory(0);
 326             fHistoQArray.AddAt(hNew,iSec);
 327         }
 328     }
 329
 330 }
 331 //_____________________________________________________________________
 332 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const  AliTPCCalibPulser &source)
 333 {
 334   //
 335   // assignment operator
 336   //
 337   if (&source == this) return *this;
 338   new (this) AliTPCCalibPulser(source);
 339
 340   return *this;
 341 }
 342 //_____________________________________________________________________
 343 AliTPCCalibPulser::~AliTPCCalibPulser()
 344 {
 345     //
 346     // destructor
 347     //
 348
 349     fCalRocArrayT0.Delete();
 350     fCalRocArrayQ.Delete();
 351     fCalRocArrayRMS.Delete();
 352     fCalRocArrayOutliers.Delete();
 353
 354     fHistoQArray.Delete();
 355     fHistoT0Array.Delete();
 356     fHistoRMSArray.Delete();
 357
 358     fPadTimesArrayEvent.Delete();
 359     fPadQArrayEvent.Delete();
 360     fPadRMSArrayEvent.Delete();
 361     fPadPedestalArrayEvent.Delete();
 362
 363     if ( fDebugStreamer) delete fDebugStreamer;
 364     delete fROC;
 365     delete fParam;
 366 }
 367 //_____________________________________________________________________
 368 Int_t AliTPCCalibPulser::Update(const Int_t icsector, /*FOLD00*/
 369                                 const Int_t icRow,
 370                                 const Int_t icPad,
 371                                 const Int_t icTimeBin,
 372                                 const Float_t csignal)
 373 {
 374     //
 375     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 376     // no extra analysis necessary. Assumes knowledge of the signal shape!
 377     // assumes that it is looped over consecutive time bins of one pad
 378     //
 379     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 380
 381     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 382
 383     //init first pad and sector in this event
 384     if ( fCurrentChannel == -1 ) {
 385         fCurrentChannel = iChannel;
 386         fCurrentSector  = icsector;
 387         fCurrentRow     = icRow;
 388     }
 389
 390     //process last pad if we change to a new one
 391     if ( iChannel != fCurrentChannel ){
 392         ProcessPad();
 393         fCurrentChannel = iChannel;
 394         fCurrentSector  = icsector;
 395         fCurrentRow     = icRow;
 396     }
 397
 398     //fill signals for current pad
 399     fPadSignal[icTimeBin]=csignal;
 400     if ( csignal > fMaxPadSignal ){
 401         fMaxPadSignal = csignal;
 402         fMaxTimeBin   = icTimeBin;
 403     }
 404     return 0;
 405 }
 406 //_____________________________________________________________________
 407 void AliTPCCalibPulser::FindPedestal(Float_t part)
 408 {
 409     //
 410     // find pedestal and noise for the current pad. Use either database or
 411     // truncated mean with part*100%
 412     //
 413     Bool_t noPedestal = kTRUE;;
 414     if (fPedestalTPC&&fPadNoiseTPC){
 415         //use pedestal database
 416         //only load new pedestals if the sector has changed
 417         if ( fCurrentSector!=fLastSector ){
 418             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 419             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 420             fLastSector=fCurrentSector;
 421         }
 422
 423         if ( fPedestalROC&&fPadNoiseROC ){
 424             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
 425             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 426             noPedestal   = kFALSE;
 427         }
 428
 429     }
 430
 431     //if we are not running with pedestal database, or for the current sector there is no information
 432     //available, calculate the pedestal and noise on the fly
 433     if ( noPedestal ) {
 434         const Int_t kPedMax = 100;  //maximum pedestal value
 435         Float_t  max    =  0;
 436         Float_t  maxPos =  0;
 437         Int_t    median =  -1;
 438         Int_t    count0 =  0;
 439         Int_t    count1 =  0;
 440         //
 441         Float_t padSignal=0;
 442         //
 443         UShort_t histo[kPedMax];
 444         memset(histo,0,kPedMax*sizeof(UShort_t));
 445
 446         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 447             padSignal = fPadSignal.GetMatrixArray()[i];
 448             if (padSignal<=0) continue;
 449             if (padSignal>max && i>10) {
 450                 max = padSignal;
 451                 maxPos = i;
 452             }
 453             if (padSignal>kPedMax-1) continue;
 454             histo[Int_t(padSignal+0.5)]++;
 455             count0++;
 456         }
 457             //
 458         for (Int_t i=1; i<kPedMax; ++i){
 459             if (count1<count0*0.5) median=i;
 460             count1+=histo[i];
 461         }
 462         // truncated mean
 463         //
 464         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 465         //
 466         for (Int_t idelta=1; idelta<10; ++idelta){
 467             if (median-idelta<=0) continue;
 468             if (median+idelta>kPedMax) continue;
 469             if (count<part*count1){
 470                 count+=histo[median-idelta];
 471                 mean +=histo[median-idelta]*(median-idelta);
 472                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 473                 count+=histo[median+idelta];
 474                 mean +=histo[median+idelta]*(median+idelta);
 475                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 476             }
 477         }
 478         fPadPedestal = 0;
 479         fPadNoise    = 0;
 480         if ( count > 0 ) {
 481             mean/=count;
 482             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 483             fPadPedestal = mean;
 484             fPadNoise    = rms;
 485         }
 486     }
 487 }
 488 //_____________________________________________________________________
 489 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
 490 {
 491     //
 492     //  Find position, signal width and height of the CE signal (last signal)
 493     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 494     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 495     //
 496
 497     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 498     Int_t   cemaxpos       = fMaxTimeBin;
 499     Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
 500     const Int_t    kCemin  = 2;             // range for the analysis of the ce signal +- channels from the peak
 501     const Int_t    kCemax  = 7;
 502
 503     if (cemaxpos!=0){
 504         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 505         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
 506             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 507             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
 508                 ceTime+=signal*(i+0.5);
 509                 ceRMS +=signal*(i+0.5)*(i+0.5);
 510                 ceQsum+=signal;
 511             }
 512         }
 513     }
 514     if (ceQmax&&ceQsum>ceSumThreshold) {
 515         ceTime/=ceQsum;
 516         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 517         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 518         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 519
 520         //Normalise Q to the pad area
 521         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
 522
 523         ceQsum/=norm;
 524     } else {
 525         ceQmax=0;
 526         ceTime=0;
 527         ceRMS =0;
 528         ceQsum=0;
 529     }
 530     param[0] = ceQmax;
 531     param[1] = ceTime;
 532     param[2] = ceRMS;
 533     qSum     = ceQsum;
 534 }
 535 //_____________________________________________________________________
 536 void AliTPCCalibPulser::ProcessPad() /*FOLD00*/
 537 {
 538     //
 539     //  Process data of current pad
 540     //
 541
 542     FindPedestal();
 543     TVectorD param(3);
 544     Float_t  Qsum;
 545     FindPulserSignal(param, Qsum);
 546
 547     Double_t meanT  = param[1];
 548     Double_t sigmaT = param[2];
 549
 550     //Fill Event T0 counter
 551     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 552
 553     //Fill Q histogram
 554     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum), fCurrentChannel );
 555
 556     //Fill RMS histogram
 557     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 558
 559
 560     //Fill debugging info
 561     if ( fDebugLevel>0 ){
 562         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
 563         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
 564         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=Qsum;
 565     }
 566
 567     ResetPad();
 568 }
 569 //_____________________________________________________________________
 570 void AliTPCCalibPulser::EndEvent() /*FOLD00*/
 571 {
 572     //
 573     //  Process data of current event
 574     //
 575
 576     //check if last pad has allready been processed, if not do so
 577     if ( fMaxTimeBin>-1 ) ProcessPad();
 578
 579     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
 580     for ( Int_t iSec = 0; iSec<72; ++iSec ){
 581         TVectorF *vTimes = GetPadTimesEvent(iSec);
 582         if ( !vTimes ) continue;
 583
 584         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
 585             Float_t Time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
 586             Float_t Time  = (*vTimes).GetMatrixArray()[iChannel];
 587
 588             GetHistoT0(iSec,kTRUE)->Fill( Time-Time0,iChannel );
 589
 590
 591             //Debug start
 592             if ( fDebugLevel>0 ){
 593                 if ( !fDebugStreamer ) {
 594                         //debug stream
 595                     TDirectory *backup = gDirectory;
 596                     fDebugStreamer = new TTreeSRedirector("deb2.root");
 597                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 598                 }
 599
 600                 Int_t row=0;
 601                 Int_t pad=0;
 602                 Int_t padc=0;
 603
 604                 Float_t Q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
 605                 Float_t RMS = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
 606
 607                 UInt_t channel=iChannel;
 608                 Int_t sector=iSec;
 609
 610                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 611                 pad = channel-fROC->GetRowIndexes(sector)[row];
 612                 padc = pad-(fROC->GetNPads(sector,row)/2);
 613
 614                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
 615                                     Form("hSignalD%d.%d.%d",sector,row,pad),
 616                                     fLastTimeBin-fFirstTimeBin,
 617                                     fFirstTimeBin,fLastTimeBin);
 618                 h1->SetDirectory(0);
 619
 620                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 621                     h1->Fill(i,fPadSignal(i));
 622
 623                 (*fDebugStreamer) << "DataPad" <<
 624                     "Event=" << fEvent <<
 625                     "Sector="<< sector <<
 626                     "Row="   << row<<
 627                     "Pad="   << pad <<
 628                     "PadC="  << padc <<
 629                     "PadSec="<< channel <<
 630                     "Time0="  << Time0 <<
 631                     "Time="  << Time <<
 632                     "RMS="   << RMS <<
 633                     "Sum="   << Q <<
 634                     "hist.=" << h1 <<
 635                     "\n";
 636
 637                 delete h1;
 638             }
 639             //Debug end
 640
 641         }
 642     }
 643
 644 }
 645 //_____________________________________________________________________
 646 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
 647 {
 648   //
 649   // Event Processing loop - AliTPCRawStream
 650   //
 651
 652   rawStream->SetOldRCUFormat(fOldRCUformat);
 653
 654   ResetEvent();
 655
 656   Bool_t withInput = kFALSE;
 657
 658   while (rawStream->Next()) {
 659
 660       Int_t isector  = rawStream->GetSector();                       //  current sector
 661       Int_t iRow     = rawStream->GetRow();                          //  current row
 662       Int_t iPad     = rawStream->GetPad();                          //  current pad
 663       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 664       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 665
 666       Update(isector,iRow,iPad,iTimeBin,signal);
 667       withInput = kTRUE;
 668   }
 669
 670   if (withInput){
 671       EndEvent();
 672   }
 673
 674   return withInput;
 675 }
 676 //_____________________________________________________________________
 677 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
 678 {
 679   //
 680   //  Event processing loop - AliRawReader
 681   //
 682
 683
 684   AliTPCRawStream rawStream(rawReader);
 685
 686   rawReader->Select("TPC");
 687
 688   return ProcessEvent(&rawStream);
 689 }
 690 //_____________________________________________________________________
 691 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
 692 {
 693   //
 694   //  Event processing loop - date event
 695   //
 696     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 697     Bool_t result=ProcessEvent(rawReader);
 698     delete rawReader;
 699     return result;
 700
 701 }
 702 //_____________________________________________________________________
 703 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
 704                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 705                                   Char_t *type, Bool_t force)
 706 {
 707     //
 708     // return pointer to Q histogram
 709     // if force is true create a new histogram if it doesn't exist allready
 710     //
 711     if ( !force || arr->UncheckedAt(sector) )
 712         return (TH2S*)arr->UncheckedAt(sector);
 713
 714     // if we are forced and histogram doesn't yes exist create it
 715     Char_t name[255], title[255];
 716
 717     sprintf(name,"hCalib%s%.2d",type,sector);
 718     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 719
 720     // new histogram with Q calib information. One value for each pad!
 721     TH2S* hist = new TH2S(name,title,
 722                           nbinsY, ymin, ymax,
 723                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 724     hist->SetDirectory(0);
 725     arr->AddAt(hist,sector);
 726     return hist;
 727 }
 728 //_____________________________________________________________________
 729 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
 730 {
 731     //
 732     // return pointer to T0 histogram
 733     // if force is true create a new histogram if it doesn't exist allready
 734     //
 735     TObjArray *arr = &fHistoT0Array;
 736     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 737 }
 738 //_____________________________________________________________________
 739 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
 740 {
 741     //
 742     // return pointer to Q histogram
 743     // if force is true create a new histogram if it doesn't exist allready
 744     //
 745     TObjArray *arr = &fHistoQArray;
 746     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 747 }
 748 //_____________________________________________________________________
 749 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
 750 {
 751     //
 752     // return pointer to Q histogram
 753     // if force is true create a new histogram if it doesn't exist allready
 754     //
 755     TObjArray *arr = &fHistoRMSArray;
 756     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 757 }
 758 //_____________________________________________________________________
 759 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
 760 {
 761     //
 762     // return pointer to Pad Info from 'arr' for the current event and sector
 763     // if force is true create it if it doesn't exist allready
 764     //
 765     if ( !force || arr->UncheckedAt(sector) )
 766         return (TVectorF*)arr->UncheckedAt(sector);
 767
 768     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 769     arr->AddAt(vect,sector);
 770     return vect;
 771 }
 772 //_____________________________________________________________________
 773 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
 774 {
 775     //
 776     // return pointer to Pad Times Array for the current event and sector
 777     // if force is true create it if it doesn't exist allready
 778     //
 779     TObjArray *arr = &fPadTimesArrayEvent;
 780     return GetPadInfoEvent(sector,arr,force);
 781 }
 782 //_____________________________________________________________________
 783 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
 784 {
 785     //
 786     // return pointer to Pad Q Array for the current event and sector
 787     // if force is true create it if it doesn't exist allready
 788     // for debugging purposes only
 789     //
 790
 791     TObjArray *arr = &fPadQArrayEvent;
 792     return GetPadInfoEvent(sector,arr,force);
 793 }
 794 //_____________________________________________________________________
 795 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
 796 {
 797     //
 798     // return pointer to Pad RMS Array for the current event and sector
 799     // if force is true create it if it doesn't exist allready
 800     // for debugging purposes only
 801     //
 802     TObjArray *arr = &fPadRMSArrayEvent;
 803     return GetPadInfoEvent(sector,arr,force);
 804 }
 805 //_____________________________________________________________________
 806 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
 807 {
 808     //
 809     // return pointer to Pad RMS Array for the current event and sector
 810     // if force is true create it if it doesn't exist allready
 811     // for debugging purposes only
 812     //
 813     TObjArray *arr = &fPadPedestalArrayEvent;
 814     return GetPadInfoEvent(sector,arr,force);
 815 }
 816 //_____________________________________________________________________
 817 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
 818 {
 819     //
 820     // return pointer to ROC Calibration
 821     // if force is true create a new histogram if it doesn't exist allready
 822     //
 823     if ( !force || arr->UncheckedAt(sector) )
 824         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 825
 826     // if we are forced and histogram doesn't yes exist create it
 827
 828     // new AliTPCCalROC for T0 information. One value for each pad!
 829     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 830     arr->AddAt(croc,sector);
 831     return croc;
 832 }
 833 //_____________________________________________________________________
 834 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
 835 {
 836     //
 837     // return pointer to Carge ROC Calibration
 838     // if force is true create a new histogram if it doesn't exist allready
 839     //
 840     TObjArray *arr = &fCalRocArrayT0;
 841     return GetCalRoc(sector, arr, force);
 842 }
 843 //_____________________________________________________________________
 844 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
 845 {
 846     //
 847     // return pointer to T0 ROC Calibration
 848     // if force is true create a new histogram if it doesn't exist allready
 849     //
 850     TObjArray *arr = &fCalRocArrayQ;
 851     return GetCalRoc(sector, arr, force);
 852 }
 853 //_____________________________________________________________________
 854 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
 855 {
 856     //
 857     // return pointer to signal width ROC Calibration
 858     // if force is true create a new histogram if it doesn't exist allready
 859     //
 860     TObjArray *arr = &fCalRocArrayRMS;
 861     return GetCalRoc(sector, arr, force);
 862 }
 863 //_____________________________________________________________________
 864 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
 865 {
 866     //
 867     // return pointer to Outliers
 868     // if force is true create a new histogram if it doesn't exist allready
 869     //
 870     TObjArray *arr = &fCalRocArrayOutliers;
 871     return GetCalRoc(sector, arr, force);
 872 }
 873 //_____________________________________________________________________
 874 void AliTPCCalibPulser::ResetEvent() /*FOLD00*/
 875 {
 876     //
 877     //  Reset global counters  -- Should be called before each event is processed
 878     //
 879     fLastSector=-1;
 880     fCurrentSector=-1;
 881     fCurrentRow=-1;
 882     fCurrentChannel=-1;
 883
 884     ResetPad();
 885
 886     fPadTimesArrayEvent.Delete();
 887     fPadQArrayEvent.Delete();
 888     fPadRMSArrayEvent.Delete();
 889     fPadPedestalArrayEvent.Delete();
 890
 891     for ( Int_t i=0; i<72; ++i ){
 892         fVTime0Offset[i]=0;
 893         fVTime0OffsetCounter[i]=0;
 894     }
 895 }
 896 //_____________________________________________________________________
 897 void AliTPCCalibPulser::ResetPad() /*FOLD00*/
 898 {
 899     //
 900     //  Reset pad infos -- Should be called after a pad has been processed
 901     //
 902     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 903         fPadSignal[i] = 0;
 904     fMaxTimeBin = -1;
 905     fMaxPadSignal = -1;
 906     fPadPedestal  = -1;
 907     fPadNoise     = -1;
 908 }
 909 //_____________________________________________________________________
 910 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
 911 {
 912     //
 913     //  Merge reference histograms of sig to the current AliTPCCalibPulser
 914     //
 915
 916     //merge histograms
 917     for (Int_t iSec=0; iSec<72; ++iSec){
 918         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
 919         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
 920         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
 921
 922
 923         if ( hRefQmerge ){
 924             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
 925             TH2S *hRefQ   = GetHistoQ(iSec);
 926             if ( hRefQ ) hRefQ->Add(hRefQmerge);
 927             else {
 928                 TH2S *hist = new TH2S(*hRefQmerge);
 929                 hist->SetDirectory(0);
 930                 fHistoQArray.AddAt(hist, iSec);
 931             }
 932             hRefQmerge->SetDirectory(dir);
 933         }
 934         if ( hRefT0merge ){
 935             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
 936             TH2S *hRefT0  = GetHistoT0(iSec);
 937             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
 938             else {
 939                 TH2S *hist = new TH2S(*hRefT0merge);
 940                 hist->SetDirectory(0);
 941                 fHistoT0Array.AddAt(hist, iSec);
 942             }
 943             hRefT0merge->SetDirectory(dir);
 944         }
 945         if ( hRefRMSmerge ){
 946             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
 947             TH2S *hRefRMS = GetHistoRMS(iSec);
 948             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
 949             else {
 950                 TH2S *hist = new TH2S(*hRefRMSmerge);
 951                 hist->SetDirectory(0);
 952                 fHistoRMSArray.AddAt(hist, iSec);
 953             }
 954             hRefRMSmerge->SetDirectory(dir);
 955         }
 956
 957     }
 958 }
 959 //_____________________________________________________________________
 960 void AliTPCCalibPulser::Analyse() /*FOLD00*/
 961 {
 962     //
 963     //  Calculate calibration constants
 964     //
 965
 966     TVectorD paramQ(3);
 967     TVectorD paramT0(3);
 968     TVectorD paramRMS(3);
 969     TMatrixD dummy(3,3);
 970
 971     for (Int_t iSec=0; iSec<72; ++iSec){
 972         TH2S *hT0 = GetHistoT0(iSec);
 973         if (!hT0 ) continue;
 974
 975         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
 976         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
 977         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
 978         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
 979
 980         TH2S *hQ   = GetHistoQ(iSec);
 981         TH2S *hRMS = GetHistoRMS(iSec);
 982
 983         Short_t *array_hQ   = hQ->GetArray();
 984         Short_t *array_hT0  = hT0->GetArray();
 985         Short_t *array_hRMS = hRMS->GetArray();
 986
 987         UInt_t nChannels = fROC->GetNChannels(iSec);
 988
 989         //debug
 990         Int_t row=0;
 991         Int_t pad=0;
 992         Int_t padc=0;
 993         //! debug
 994
 995         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
 996
 997
 998             Float_t cogTime0 = -1000;
 999             Float_t cogQ     = -1000;
1000             Float_t cogRMS   = -1000;
1001             Float_t cogOut   = 0;
1002
1003
1004             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1005             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1006             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1007
1008 /*
1009             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1010             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1011             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1012             cogQ     = paramQ[1];
1013             cogTime0 = paramT0[1];
1014             cogRMS   = paramRMS[1];
1015 */
1016             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1017             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1018             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1019
1020
1021
1022             /*
1023             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1024                 cogOut = 1;
1025                 cogTime0 = 0;
1026                 cogQ     = 0;
1027                 cogRMS   = 0;
1028             }
1029 */
1030             rocQ->SetValue(iChannel, cogQ*cogQ);
1031             rocT0->SetValue(iChannel, cogTime0);
1032             rocRMS->SetValue(iChannel, cogRMS);
1033             rocOut->SetValue(iChannel, cogOut);
1034
1035
1036             //debug
1037             if ( fDebugLevel > 0 ){
1038                 if ( !fDebugStreamer ) {
1039                         //debug stream
1040                     TDirectory *backup = gDirectory;
1041                     fDebugStreamer = new TTreeSRedirector("deb2.root");
1042                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1043                 }
1044
1045                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1046                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1047                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1048
1049                 (*fDebugStreamer) << "DataEnd" <<
1050                     "Sector="  << iSec      <<
1051                     "Pad="     << pad       <<
1052                     "PadC="    << padc      <<
1053                     "Row="     << row       <<
1054                     "PadSec="  << iChannel   <<
1055                     "Q="       << cogQ      <<
1056                     "T0="      << cogTime0  <<
1057                     "RMS="     << cogRMS    <<
1058                     "\n";
1059             }
1060             //! debug
1061
1062         }
1063
1064     }
1065     delete fDebugStreamer;
1066     fDebugStreamer = 0x0;
1067 }
1068 //_____________________________________________________________________
1069 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1070 {
1071     //
1072     //  Write class to file
1073     //
1074
1075     TString sDir(dir);
1076     TString option;
1077
1078     if ( append )
1079         option = "update";
1080     else
1081         option = "recreate";
1082
1083     TDirectory *backup = gDirectory;
1084     TFile f(filename,option.Data());
1085     f.cd();
1086     if ( !sDir.IsNull() ){
1087         f.mkdir(sDir.Data());
1088         f.cd(sDir);
1089     }
1090     this->Write();
1091     f.Close();
1092
1093     if ( backup ) backup->cd();
1094 }
1095 //_____________________________________________________________________
1096 //_________________________  Test Functions ___________________________
1097 //_____________________________________________________________________
1098 TObjArray* AliTPCCalibPulser::TestBinning()
1099 {
1100     //
1101     //  Function to test the binning of the reference histograms
1102     //  type: T0, Q or RMS
1103     //  mode: 0 - number of filled bins per channel
1104     //        1 - number of empty bins between filled bins in one ROC
1105     //  returns TObjArray with the test histograms type*2+mode:
1106     //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1107
1108
1109     TObjArray *histArray = new TObjArray(6);
1110     const Char_t *type[] = {"T0","Q","RMS"};
1111     Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1112
1113     for (Int_t itype = 0; itype<3; ++itype){
1114         for (Int_t imode=0; imode<2; ++imode){
1115             Int_t icount = itype*2+imode;
1116             histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1117                                       Form("Test Binning of '%s', mode - %d",type[itype],imode),
1118                                       72,0,72),
1119                              icount);
1120         }
1121     }
1122
1123
1124     TH2S *hRef=0x0;
1125     Short_t *array=0x0;
1126     for (Int_t itype = 0; itype<3; ++itype){
1127         for (Int_t iSec=0; iSec<72; ++iSec){
1128             if ( itype == 0 ) hRef = GetHistoT0(iSec);
1129             if ( itype == 1 ) hRef = GetHistoQ(iSec);
1130             if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1131             if ( hRef == 0x0 ) continue;
1132             array = (hRef->GetArray());
1133             UInt_t nChannels = fROC->GetNChannels(iSec);
1134
1135             Int_t nempty=0;
1136             for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1137                 Int_t nfilled=0;
1138                 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1139                 Int_t c1 = 0;
1140                 Int_t c2 = 0;
1141                 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1142                     if ( array[offset+iBin]>0 ) {
1143                         nfilled++;
1144                         if ( c1 && c2 ) nempty++;
1145                         else c1 = 1;
1146                     }
1147                     else if ( c1 ) c2 = 1;
1148
1149                 }
1150                 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1151             }
1152             ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1153         }
1154     }
1155     return histArray;
1156 }