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