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 fParam;
 365 }
 366 //_____________________________________________________________________
 367 Int_t AliTPCCalibPulser::Update(const Int_t icsector, /*FOLD00*/
 368                                 const Int_t icRow,
 369                                 const Int_t icPad,
 370                                 const Int_t icTimeBin,
 371                                 const Float_t csignal)
 372 {
 373     //
 374     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 375     // no extra analysis necessary. Assumes knowledge of the signal shape!
 376     // assumes that it is looped over consecutive time bins of one pad
 377     //
 378     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 379
 380     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 381
 382     //init first pad and sector in this event
 383     if ( fCurrentChannel == -1 ) {
 384         fCurrentChannel = iChannel;
 385         fCurrentSector  = icsector;
 386         fCurrentRow     = icRow;
 387     }
 388
 389     //process last pad if we change to a new one
 390     if ( iChannel != fCurrentChannel ){
 391         ProcessPad();
 392         fCurrentChannel = iChannel;
 393         fCurrentSector  = icsector;
 394         fCurrentRow     = icRow;
 395     }
 396
 397     //fill signals for current pad
 398     fPadSignal[icTimeBin]=csignal;
 399     if ( csignal > fMaxPadSignal ){
 400         fMaxPadSignal = csignal;
 401         fMaxTimeBin   = icTimeBin;
 402     }
 403     return 0;
 404 }
 405 //_____________________________________________________________________
 406 void AliTPCCalibPulser::FindPedestal(Float_t part)
 407 {
 408     //
 409     // find pedestal and noise for the current pad. Use either database or
 410     // truncated mean with part*100%
 411     //
 412     Bool_t noPedestal = kTRUE;;
 413     if (fPedestalTPC&&fPadNoiseTPC){
 414         //use pedestal database
 415         //only load new pedestals if the sector has changed
 416         if ( fCurrentSector!=fLastSector ){
 417             fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 418             fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
 419             fLastSector=fCurrentSector;
 420         }
 421
 422         if ( fPedestalROC&&fPadNoiseROC ){
 423             fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
 424             fPadNoise    = fPadNoiseROC->GetValue(fCurrentChannel);
 425             noPedestal   = kFALSE;
 426         }
 427
 428     }
 429
 430     //if we are not running with pedestal database, or for the current sector there is no information
 431     //available, calculate the pedestal and noise on the fly
 432     if ( noPedestal ) {
 433         const Int_t kPedMax = 100;  //maximum pedestal value
 434         Float_t  max    =  0;
 435         Float_t  maxPos =  0;
 436         Int_t    median =  -1;
 437         Int_t    count0 =  0;
 438         Int_t    count1 =  0;
 439         //
 440         Float_t padSignal=0;
 441         //
 442         UShort_t histo[kPedMax];
 443         memset(histo,0,kPedMax*sizeof(UShort_t));
 444
 445         for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
 446             padSignal = fPadSignal.GetMatrixArray()[i];
 447             if (padSignal<=0) continue;
 448             if (padSignal>max && i>10) {
 449                 max = padSignal;
 450                 maxPos = i;
 451             }
 452             if (padSignal>kPedMax-1) continue;
 453             histo[Int_t(padSignal+0.5)]++;
 454             count0++;
 455         }
 456             //
 457         for (Int_t i=1; i<kPedMax; ++i){
 458             if (count1<count0*0.5) median=i;
 459             count1+=histo[i];
 460         }
 461         // truncated mean
 462         //
 463         Float_t count=histo[median] ,mean=histo[median]*median,  rms=histo[median]*median*median ;
 464         //
 465         for (Int_t idelta=1; idelta<10; ++idelta){
 466             if (median-idelta<=0) continue;
 467             if (median+idelta>kPedMax) continue;
 468             if (count<part*count1){
 469                 count+=histo[median-idelta];
 470                 mean +=histo[median-idelta]*(median-idelta);
 471                 rms  +=histo[median-idelta]*(median-idelta)*(median-idelta);
 472                 count+=histo[median+idelta];
 473                 mean +=histo[median+idelta]*(median+idelta);
 474                 rms  +=histo[median+idelta]*(median+idelta)*(median+idelta);
 475             }
 476         }
 477         fPadPedestal = 0;
 478         fPadNoise    = 0;
 479         if ( count > 0 ) {
 480             mean/=count;
 481             rms    = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
 482             fPadPedestal = mean;
 483             fPadNoise    = rms;
 484         }
 485     }
 486 }
 487 //_____________________________________________________________________
 488 void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
 489 {
 490     //
 491     //  Find position, signal width and height of the CE signal (last signal)
 492     //  param[0] = Qmax, param[1] = mean time, param[2] = rms;
 493     //  maxima: array of local maxima of the pad signal use the one closest to the mean CE position
 494     //
 495
 496     Float_t ceQmax  =0, ceQsum=0, ceTime=0, ceRMS=0;
 497     Int_t   cemaxpos       = fMaxTimeBin;
 498     Float_t ceSumThreshold = 8.*fPadNoise;  // threshold for the signal sum
 499     const Int_t    kCemin  = 2;             // range for the analysis of the ce signal +- channels from the peak
 500     const Int_t    kCemax  = 7;
 501
 502     if (cemaxpos!=0){
 503         ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
 504         for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
 505             Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
 506             if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
 507                 ceTime+=signal*(i+0.5);
 508                 ceRMS +=signal*(i+0.5)*(i+0.5);
 509                 ceQsum+=signal;
 510             }
 511         }
 512     }
 513     if (ceQmax&&ceQsum>ceSumThreshold) {
 514         ceTime/=ceQsum;
 515         ceRMS  = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
 516         fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime;   // mean time for each sector
 517         fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
 518
 519         //Normalise Q to the pad area
 520         Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
 521
 522         ceQsum/=norm;
 523     } else {
 524         ceQmax=0;
 525         ceTime=0;
 526         ceRMS =0;
 527         ceQsum=0;
 528     }
 529     param[0] = ceQmax;
 530     param[1] = ceTime;
 531     param[2] = ceRMS;
 532     qSum     = ceQsum;
 533 }
 534 //_____________________________________________________________________
 535 void AliTPCCalibPulser::ProcessPad() /*FOLD00*/
 536 {
 537     //
 538     //  Process data of current pad
 539     //
 540
 541     FindPedestal();
 542     TVectorD param(3);
 543     Float_t  Qsum;
 544     FindPulserSignal(param, Qsum);
 545
 546     Double_t meanT  = param[1];
 547     Double_t sigmaT = param[2];
 548
 549     //Fill Event T0 counter
 550     (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
 551
 552     //Fill Q histogram
 553     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum), fCurrentChannel );
 554
 555     //Fill RMS histogram
 556     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 557
 558
 559     //Fill debugging info
 560     if ( fDebugLevel>0 ){
 561         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
 562         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
 563         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=Qsum;
 564     }
 565
 566     ResetPad();
 567 }
 568 //_____________________________________________________________________
 569 void AliTPCCalibPulser::EndEvent() /*FOLD00*/
 570 {
 571     //
 572     //  Process data of current event
 573     //
 574
 575     //check if last pad has allready been processed, if not do so
 576     if ( fMaxTimeBin>-1 ) ProcessPad();
 577
 578     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
 579     for ( Int_t iSec = 0; iSec<72; ++iSec ){
 580         TVectorF *vTimes = GetPadTimesEvent(iSec);
 581         if ( !vTimes ) continue;
 582
 583         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
 584             Float_t Time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
 585             Float_t Time  = (*vTimes).GetMatrixArray()[iChannel];
 586
 587             GetHistoT0(iSec,kTRUE)->Fill( Time-Time0,iChannel );
 588
 589
 590             //Debug start
 591             if ( fDebugLevel>0 ){
 592                 if ( !fDebugStreamer ) {
 593                         //debug stream
 594                     TDirectory *backup = gDirectory;
 595                     fDebugStreamer = new TTreeSRedirector("deb2.root");
 596                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 597                 }
 598
 599                 Int_t row=0;
 600                 Int_t pad=0;
 601                 Int_t padc=0;
 602
 603                 Float_t Q   = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
 604                 Float_t RMS = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
 605
 606                 UInt_t channel=iChannel;
 607                 Int_t sector=iSec;
 608
 609                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 610                 pad = channel-fROC->GetRowIndexes(sector)[row];
 611                 padc = pad-(fROC->GetNPads(sector,row)/2);
 612
 613                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
 614                                     Form("hSignalD%d.%d.%d",sector,row,pad),
 615                                     fLastTimeBin-fFirstTimeBin,
 616                                     fFirstTimeBin,fLastTimeBin);
 617                 h1->SetDirectory(0);
 618
 619                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 620                     h1->Fill(i,fPadSignal(i));
 621
 622                 (*fDebugStreamer) << "DataPad" <<
 623                     "Event=" << fEvent <<
 624                     "Sector="<< sector <<
 625                     "Row="   << row<<
 626                     "Pad="   << pad <<
 627                     "PadC="  << padc <<
 628                     "PadSec="<< channel <<
 629                     "Time0="  << Time0 <<
 630                     "Time="  << Time <<
 631                     "RMS="   << RMS <<
 632                     "Sum="   << Q <<
 633                     "hist.=" << h1 <<
 634                     "\n";
 635
 636                 delete h1;
 637             }
 638             //Debug end
 639
 640         }
 641     }
 642
 643 }
 644 //_____________________________________________________________________
 645 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
 646 {
 647   //
 648   // Event Processing loop - AliTPCRawStream
 649   //
 650
 651   rawStream->SetOldRCUFormat(fOldRCUformat);
 652
 653   ResetEvent();
 654
 655   Bool_t withInput = kFALSE;
 656
 657   while (rawStream->Next()) {
 658
 659       Int_t isector  = rawStream->GetSector();                       //  current sector
 660       Int_t iRow     = rawStream->GetRow();                          //  current row
 661       Int_t iPad     = rawStream->GetPad();                          //  current pad
 662       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 663       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 664
 665       Update(isector,iRow,iPad,iTimeBin,signal);
 666       withInput = kTRUE;
 667   }
 668
 669   if (withInput){
 670       EndEvent();
 671   }
 672
 673   return withInput;
 674 }
 675 //_____________________________________________________________________
 676 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
 677 {
 678   //
 679   //  Event processing loop - AliRawReader
 680   //
 681
 682
 683   AliTPCRawStream rawStream(rawReader);
 684
 685   rawReader->Select("TPC");
 686
 687   return ProcessEvent(&rawStream);
 688 }
 689 //_____________________________________________________________________
 690 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
 691 {
 692   //
 693   //  Event processing loop - date event
 694   //
 695     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 696     Bool_t result=ProcessEvent(rawReader);
 697     delete rawReader;
 698     return result;
 699
 700 }
 701 //_____________________________________________________________________
 702 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
 703                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 704                                   Char_t *type, Bool_t force)
 705 {
 706     //
 707     // return pointer to Q histogram
 708     // if force is true create a new histogram if it doesn't exist allready
 709     //
 710     if ( !force || arr->UncheckedAt(sector) )
 711         return (TH2S*)arr->UncheckedAt(sector);
 712
 713     // if we are forced and histogram doesn't yes exist create it
 714     Char_t name[255], title[255];
 715
 716     sprintf(name,"hCalib%s%.2d",type,sector);
 717     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 718
 719     // new histogram with Q calib information. One value for each pad!
 720     TH2S* hist = new TH2S(name,title,
 721                           nbinsY, ymin, ymax,
 722                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 723     hist->SetDirectory(0);
 724     arr->AddAt(hist,sector);
 725     return hist;
 726 }
 727 //_____________________________________________________________________
 728 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
 729 {
 730     //
 731     // return pointer to T0 histogram
 732     // if force is true create a new histogram if it doesn't exist allready
 733     //
 734     TObjArray *arr = &fHistoT0Array;
 735     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 736 }
 737 //_____________________________________________________________________
 738 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
 739 {
 740     //
 741     // return pointer to Q histogram
 742     // if force is true create a new histogram if it doesn't exist allready
 743     //
 744     TObjArray *arr = &fHistoQArray;
 745     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 746 }
 747 //_____________________________________________________________________
 748 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
 749 {
 750     //
 751     // return pointer to Q histogram
 752     // if force is true create a new histogram if it doesn't exist allready
 753     //
 754     TObjArray *arr = &fHistoRMSArray;
 755     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 756 }
 757 //_____________________________________________________________________
 758 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
 759 {
 760     //
 761     // return pointer to Pad Info from 'arr' for the current event and sector
 762     // if force is true create it if it doesn't exist allready
 763     //
 764     if ( !force || arr->UncheckedAt(sector) )
 765         return (TVectorF*)arr->UncheckedAt(sector);
 766
 767     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 768     arr->AddAt(vect,sector);
 769     return vect;
 770 }
 771 //_____________________________________________________________________
 772 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
 773 {
 774     //
 775     // return pointer to Pad Times Array for the current event and sector
 776     // if force is true create it if it doesn't exist allready
 777     //
 778     TObjArray *arr = &fPadTimesArrayEvent;
 779     return GetPadInfoEvent(sector,arr,force);
 780 }
 781 //_____________________________________________________________________
 782 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
 783 {
 784     //
 785     // return pointer to Pad Q Array for the current event and sector
 786     // if force is true create it if it doesn't exist allready
 787     // for debugging purposes only
 788     //
 789
 790     TObjArray *arr = &fPadQArrayEvent;
 791     return GetPadInfoEvent(sector,arr,force);
 792 }
 793 //_____________________________________________________________________
 794 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
 795 {
 796     //
 797     // return pointer to Pad RMS Array for the current event and sector
 798     // if force is true create it if it doesn't exist allready
 799     // for debugging purposes only
 800     //
 801     TObjArray *arr = &fPadRMSArrayEvent;
 802     return GetPadInfoEvent(sector,arr,force);
 803 }
 804 //_____________________________________________________________________
 805 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
 806 {
 807     //
 808     // return pointer to Pad RMS Array for the current event and sector
 809     // if force is true create it if it doesn't exist allready
 810     // for debugging purposes only
 811     //
 812     TObjArray *arr = &fPadPedestalArrayEvent;
 813     return GetPadInfoEvent(sector,arr,force);
 814 }
 815 //_____________________________________________________________________
 816 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
 817 {
 818     //
 819     // return pointer to ROC Calibration
 820     // if force is true create a new histogram if it doesn't exist allready
 821     //
 822     if ( !force || arr->UncheckedAt(sector) )
 823         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 824
 825     // if we are forced and histogram doesn't yes exist create it
 826
 827     // new AliTPCCalROC for T0 information. One value for each pad!
 828     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 829     arr->AddAt(croc,sector);
 830     return croc;
 831 }
 832 //_____________________________________________________________________
 833 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
 834 {
 835     //
 836     // return pointer to Carge ROC Calibration
 837     // if force is true create a new histogram if it doesn't exist allready
 838     //
 839     TObjArray *arr = &fCalRocArrayT0;
 840     return GetCalRoc(sector, arr, force);
 841 }
 842 //_____________________________________________________________________
 843 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
 844 {
 845     //
 846     // return pointer to T0 ROC Calibration
 847     // if force is true create a new histogram if it doesn't exist allready
 848     //
 849     TObjArray *arr = &fCalRocArrayQ;
 850     return GetCalRoc(sector, arr, force);
 851 }
 852 //_____________________________________________________________________
 853 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
 854 {
 855     //
 856     // return pointer to signal width ROC Calibration
 857     // if force is true create a new histogram if it doesn't exist allready
 858     //
 859     TObjArray *arr = &fCalRocArrayRMS;
 860     return GetCalRoc(sector, arr, force);
 861 }
 862 //_____________________________________________________________________
 863 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
 864 {
 865     //
 866     // return pointer to Outliers
 867     // if force is true create a new histogram if it doesn't exist allready
 868     //
 869     TObjArray *arr = &fCalRocArrayOutliers;
 870     return GetCalRoc(sector, arr, force);
 871 }
 872 //_____________________________________________________________________
 873 void AliTPCCalibPulser::ResetEvent() /*FOLD00*/
 874 {
 875     //
 876     //  Reset global counters  -- Should be called before each event is processed
 877     //
 878     fLastSector=-1;
 879     fCurrentSector=-1;
 880     fCurrentRow=-1;
 881     fCurrentChannel=-1;
 882
 883     ResetPad();
 884
 885     fPadTimesArrayEvent.Delete();
 886     fPadQArrayEvent.Delete();
 887     fPadRMSArrayEvent.Delete();
 888     fPadPedestalArrayEvent.Delete();
 889
 890     for ( Int_t i=0; i<72; ++i ){
 891         fVTime0Offset[i]=0;
 892         fVTime0OffsetCounter[i]=0;
 893     }
 894 }
 895 //_____________________________________________________________________
 896 void AliTPCCalibPulser::ResetPad() /*FOLD00*/
 897 {
 898     //
 899     //  Reset pad infos -- Should be called after a pad has been processed
 900     //
 901     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
 902         fPadSignal[i] = 0;
 903     fMaxTimeBin = -1;
 904     fMaxPadSignal = -1;
 905     fPadPedestal  = -1;
 906     fPadNoise     = -1;
 907 }
 908 //_____________________________________________________________________
 909 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
 910 {
 911     //
 912     //  Merge reference histograms of sig to the current AliTPCCalibPulser
 913     //
 914
 915     //merge histograms
 916     for (Int_t iSec=0; iSec<72; ++iSec){
 917         TH2S *hRefQmerge   = sig->GetHistoQ(iSec);
 918         TH2S *hRefT0merge  = sig->GetHistoT0(iSec);
 919         TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
 920
 921
 922         if ( hRefQmerge ){
 923             TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
 924             TH2S *hRefQ   = GetHistoQ(iSec);
 925             if ( hRefQ ) hRefQ->Add(hRefQmerge);
 926             else {
 927                 TH2S *hist = new TH2S(*hRefQmerge);
 928                 hist->SetDirectory(0);
 929                 fHistoQArray.AddAt(hist, iSec);
 930             }
 931             hRefQmerge->SetDirectory(dir);
 932         }
 933         if ( hRefT0merge ){
 934             TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
 935             TH2S *hRefT0  = GetHistoT0(iSec);
 936             if ( hRefT0 ) hRefT0->Add(hRefT0merge);
 937             else {
 938                 TH2S *hist = new TH2S(*hRefT0merge);
 939                 hist->SetDirectory(0);
 940                 fHistoT0Array.AddAt(hist, iSec);
 941             }
 942             hRefT0merge->SetDirectory(dir);
 943         }
 944         if ( hRefRMSmerge ){
 945             TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
 946             TH2S *hRefRMS = GetHistoRMS(iSec);
 947             if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
 948             else {
 949                 TH2S *hist = new TH2S(*hRefRMSmerge);
 950                 hist->SetDirectory(0);
 951                 fHistoRMSArray.AddAt(hist, iSec);
 952             }
 953             hRefRMSmerge->SetDirectory(dir);
 954         }
 955
 956     }
 957 }
 958 //_____________________________________________________________________
 959 void AliTPCCalibPulser::Analyse() /*FOLD00*/
 960 {
 961     //
 962     //  Calculate calibration constants
 963     //
 964
 965     TVectorD paramQ(3);
 966     TVectorD paramT0(3);
 967     TVectorD paramRMS(3);
 968     TMatrixD dummy(3,3);
 969
 970     for (Int_t iSec=0; iSec<72; ++iSec){
 971         TH2S *hT0 = GetHistoT0(iSec);
 972         if (!hT0 ) continue;
 973
 974         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
 975         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
 976         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
 977         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
 978
 979         TH2S *hQ   = GetHistoQ(iSec);
 980         TH2S *hRMS = GetHistoRMS(iSec);
 981
 982         Short_t *array_hQ   = hQ->GetArray();
 983         Short_t *array_hT0  = hT0->GetArray();
 984         Short_t *array_hRMS = hRMS->GetArray();
 985
 986         UInt_t nChannels = fROC->GetNChannels(iSec);
 987
 988         //debug
 989         Int_t row=0;
 990         Int_t pad=0;
 991         Int_t padc=0;
 992         //! debug
 993
 994         for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
 995
 996
 997             Float_t cogTime0 = -1000;
 998             Float_t cogQ     = -1000;
 999             Float_t cogRMS   = -1000;
1000             Float_t cogOut   = 0;
1001
1002
1003             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1004             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1005             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1006
1007 /*
1008             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1009             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1010             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1011             cogQ     = paramQ[1];
1012             cogTime0 = paramT0[1];
1013             cogRMS   = paramRMS[1];
1014 */
1015             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1016             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1017             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1018
1019
1020
1021             /*
1022             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1023                 cogOut = 1;
1024                 cogTime0 = 0;
1025                 cogQ     = 0;
1026                 cogRMS   = 0;
1027             }
1028 */
1029             rocQ->SetValue(iChannel, cogQ*cogQ);
1030             rocT0->SetValue(iChannel, cogTime0);
1031             rocRMS->SetValue(iChannel, cogRMS);
1032             rocOut->SetValue(iChannel, cogOut);
1033
1034
1035             //debug
1036             if ( fDebugLevel > 0 ){
1037                 if ( !fDebugStreamer ) {
1038                         //debug stream
1039                     TDirectory *backup = gDirectory;
1040                     fDebugStreamer = new TTreeSRedirector("deb2.root");
1041                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
1042                 }
1043
1044                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1045                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1046                 padc = pad-(fROC->GetNPads(iSec,row)/2);
1047
1048                 (*fDebugStreamer) << "DataEnd" <<
1049                     "Sector="  << iSec      <<
1050                     "Pad="     << pad       <<
1051                     "PadC="    << padc      <<
1052                     "Row="     << row       <<
1053                     "PadSec="  << iChannel   <<
1054                     "Q="       << cogQ      <<
1055                     "T0="      << cogTime0  <<
1056                     "RMS="     << cogRMS    <<
1057                     "\n";
1058             }
1059             //! debug
1060
1061         }
1062
1063     }
1064     delete fDebugStreamer;
1065     fDebugStreamer = 0x0;
1066 }
1067 //_____________________________________________________________________
1068 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1069 {
1070     //
1071     //  Write class to file
1072     //
1073
1074     TString sDir(dir);
1075     TString option;
1076
1077     if ( append )
1078         option = "update";
1079     else
1080         option = "recreate";
1081
1082     TDirectory *backup = gDirectory;
1083     TFile f(filename,option.Data());
1084     f.cd();
1085     if ( !sDir.IsNull() ){
1086         f.mkdir(sDir.Data());
1087         f.cd(sDir);
1088     }
1089     this->Write();
1090     f.Close();
1091
1092     if ( backup ) backup->cd();
1093 }
1094 //_____________________________________________________________________
1095 //_________________________  Test Functions ___________________________
1096 //_____________________________________________________________________
1097 TObjArray* AliTPCCalibPulser::TestBinning()
1098 {
1099     //
1100     //  Function to test the binning of the reference histograms
1101     //  type: T0, Q or RMS
1102     //  mode: 0 - number of filled bins per channel
1103     //        1 - number of empty bins between filled bins in one ROC
1104     //  returns TObjArray with the test histograms type*2+mode:
1105     //  position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1106
1107
1108     TObjArray *histArray = new TObjArray(6);
1109     const Char_t *type[] = {"T0","Q","RMS"};
1110     Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1111
1112     for (Int_t itype = 0; itype<3; ++itype){
1113         for (Int_t imode=0; imode<2; ++imode){
1114             Int_t icount = itype*2+imode;
1115             histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1116                                       Form("Test Binning of '%s', mode - %d",type[itype],imode),
1117                                       72,0,72),
1118                              icount);
1119         }
1120     }
1121
1122
1123     TH2S *hRef=0x0;
1124     Short_t *array=0x0;
1125     for (Int_t itype = 0; itype<3; ++itype){
1126         for (Int_t iSec=0; iSec<72; ++iSec){
1127             if ( itype == 0 ) hRef = GetHistoT0(iSec);
1128             if ( itype == 1 ) hRef = GetHistoQ(iSec);
1129             if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1130             if ( hRef == 0x0 ) continue;
1131             array = (hRef->GetArray());
1132             UInt_t nChannels = fROC->GetNChannels(iSec);
1133
1134             Int_t nempty=0;
1135             for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1136                 Int_t nfilled=0;
1137                 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1138                 Int_t c1 = 0;
1139                 Int_t c2 = 0;
1140                 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1141                     if ( array[offset+iBin]>0 ) {
1142                         nfilled++;
1143                         if ( c1 && c2 ) nempty++;
1144                         else c1 = 1;
1145                     }
1146                     else if ( c1 ) c2 = 1;
1147
1148                 }
1149                 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1150             }
1151             ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1152         }
1153     }
1154     return histArray;
1155 }