TPC/AliTPCCalibSignal.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17
  18
  19
  20 //-------------------------------------------------------
  21 //          Implementation of the TPC pulser calibration
  22 //
  23 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
  24 //
  25 //
  26 //-------------------------------------------------------
  27
  28
  29 /* $Id$ */
  30
  31
  32
  33 //Root includes
  34 #include <TObjArray.h>
  35 #include <TH1F.h>
  36 #include <TH2S.h>
  37 #include <TString.h>
  38 #include <TVectorF.h>
  39 #include <TMath.h>
  40
  41 #include <TDirectory.h>
  42 #include <TSystem.h>
  43 #include <TFile.h>
  44
  45 //AliRoot includes
  46 #include "AliRawReader.h"
  47 #include "AliRawReaderRoot.h"
  48 #include "AliRawReaderDate.h"
  49 #include "AliTPCRawStream.h"
  50 #include "AliTPCCalROC.h"
  51 #include "AliTPCCalPad.h"
  52 #include "AliTPCROC.h"
  53 #include "AliTPCParam.h"
  54 #include "AliTPCCalibSignal.h"
  55 #include "AliTPCcalibDB.h"
  56 #include "AliMathBase.h"
  57 #include "TTreeStream.h"
  58
  59 //date
  60 #include "event.h"
  61 ClassImp(AliTPCCalibSignal) /*FOLD00*/
  62
  63 AliTPCCalibSignal::AliTPCCalibSignal() : /*FOLD00*/
  64     TObject(),
  65     fFirstTimeBin(60),
  66     fLastTimeBin(120),
  67     fFirstTimeBinT0(-15),
  68     fLastTimeBinT0(15),
  69     fNbinsT0(200),
  70     fXminT0(-2),
  71     fXmaxT0(2),
  72     fNbinsQ(200),
  73     fXminQ(14),
  74     fXmaxQ(55),
  75     fNbinsRMS(100),
  76     fXminRMS(0),
  77     fXmaxRMS(5),
  78     fLastSector(-1),
  79     fROC(AliTPCROC::Instance()),
  80     fParam(new AliTPCParam),
  81     fPedestalTPC(0x0),
  82     fBpedestal(kFALSE),
  83     fCalRocArrayT0(72),
  84     fCalRocArrayQ(72),
  85     fCalRocArrayRMS(72),
  86     fCalRocArrayOutliers(72),
  87     fHistoQArray(72),
  88     fHistoT0Array(72),
  89     fHistoRMSArray(72),
  90     fPadTimesArrayEvent(72),
  91     fPadQArrayEvent(72),
  92     fPadRMSArrayEvent(72),
  93     fPadPedestalArrayEvent(72),
  94     fCurrentChannel(-1),
  95     fCurrentSector(-1),
  96     fCurrentRow(-1),
  97     fMaxPadSignal(-1),
  98     fMaxTimeBin(-1),
  99     fPadSignal(1024),
 100     fVTime0Offset1(72),
 101     fVTime0Offset1Counter(72),
 102     fEvent(-1),
 103     fDebugStreamer(0x0),
 104     fDebugLevel(0)
 105 {
 106     //
 107     // AliTPCSignal default constructor
 108     //
 109
 110 }
 111 //_____________________________________________________________________
 112 AliTPCCalibSignal::AliTPCCalibSignal(const AliTPCCalibSignal &sig) :
 113     TObject(sig),
 114     fFirstTimeBin(sig.fFirstTimeBin),
 115     fLastTimeBin(sig.fLastTimeBin),
 116     fFirstTimeBinT0(sig.fFirstTimeBinT0),
 117     fLastTimeBinT0(sig.fLastTimeBinT0),
 118     fNbinsT0(sig.fNbinsT0),
 119     fXminT0(sig.fXminT0),
 120     fXmaxT0(sig.fXmaxT0),
 121     fNbinsQ(sig.fNbinsQ),
 122     fXminQ(sig.fXminQ),
 123     fXmaxQ(sig.fXmaxQ),
 124     fNbinsRMS(sig.fNbinsRMS),
 125     fXminRMS(sig.fXminRMS),
 126     fXmaxRMS(sig.fXmaxRMS),
 127     fLastSector(-1),
 128     fROC(AliTPCROC::Instance()),
 129     fParam(new AliTPCParam),
 130     fPedestalTPC(sig.fPedestalTPC),
 131     fBpedestal(sig.fBpedestal),
 132     fCalRocArrayT0(72),
 133     fCalRocArrayQ(72),
 134     fCalRocArrayRMS(72),
 135     fCalRocArrayOutliers(72),
 136     fHistoQArray(72),
 137     fHistoT0Array(72),
 138     fHistoRMSArray(72),
 139     fPadTimesArrayEvent(72),
 140     fPadQArrayEvent(72),
 141     fPadRMSArrayEvent(72),
 142     fPadPedestalArrayEvent(72),
 143     fCurrentChannel(-1),
 144     fCurrentSector(-1),
 145     fCurrentRow(-1),
 146     fMaxPadSignal(-1),
 147     fMaxTimeBin(-1),
 148     fPadSignal(1024),
 149     fVTime0Offset1(72),
 150     fVTime0Offset1Counter(72),
 151     fEvent(-1),
 152     fDebugStreamer(0x0),
 153     fDebugLevel(sig.fDebugLevel)
 154 {
 155     //
 156     // AliTPCSignal default constructor
 157     //
 158
 159     for (Int_t iSec = 0; iSec < 72; iSec++){
 160         const AliTPCCalROC *calQ   = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
 161         const AliTPCCalROC *calT0  = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
 162         const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
 163         const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
 164
 165         const TH2S *hQ   = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
 166         const TH2S *hT0  = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
 167         const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
 168
 169         if ( calQ   != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
 170         if ( calT0  != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
 171         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 172         if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
 173
 174         if ( hQ != 0x0 ){
 175             TH2S *hNew = new TH2S(*hQ);
 176             hNew->SetDirectory(0);
 177             fHistoQArray.AddAt(hNew,iSec);
 178         }
 179         if ( hT0 != 0x0 ){
 180             TH2S *hNew = new TH2S(*hT0);
 181             hNew->SetDirectory(0);
 182             fHistoQArray.AddAt(hNew,iSec);
 183         }
 184         if ( hRMS != 0x0 ){
 185             TH2S *hNew = new TH2S(*hRMS);
 186             hNew->SetDirectory(0);
 187             fHistoQArray.AddAt(hNew,iSec);
 188         }
 189     }
 190
 191 }
 192 //_____________________________________________________________________
 193 AliTPCCalibSignal& AliTPCCalibSignal::operator = (const  AliTPCCalibSignal &source)
 194 {
 195   //
 196   // assignment operator
 197   //
 198   if (&source == this) return *this;
 199   new (this) AliTPCCalibSignal(source);
 200
 201   return *this;
 202 }
 203 //_____________________________________________________________________
 204 AliTPCCalibSignal::~AliTPCCalibSignal()
 205 {
 206     //
 207     // destructor
 208     //
 209
 210     fCalRocArrayT0.Delete();
 211     fCalRocArrayQ.Delete();
 212     fCalRocArrayRMS.Delete();
 213
 214     fHistoQArray.Delete();
 215     fHistoT0Array.Delete();
 216     fHistoRMSArray.Delete();
 217
 218     fPadTimesArrayEvent.Delete();
 219     fPadQArrayEvent.Delete();
 220     fPadRMSArrayEvent.Delete();
 221     fPadPedestalArrayEvent.Delete();
 222
 223     if ( fDebugStreamer) delete fDebugStreamer;
 224     delete fROC;
 225     delete fParam;
 226 }
 227 //_____________________________________________________________________
 228 Int_t AliTPCCalibSignal::Update(const Int_t icsector, /*FOLD00*/
 229                                 const Int_t icRow,
 230                                 const Int_t icPad,
 231                                 const Int_t icTimeBin,
 232                                 const Float_t csignal)
 233 {
 234     //
 235     // Signal filling methode on the fly pedestal and Time offset correction if necessary.
 236     // no extra analysis necessary. Assumes knowledge of the signal shape!
 237     // assumes that it is looped over consecutive time bins of one pad
 238     //
 239     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 240
 241     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 242
 243     //init first pad and sector in this event
 244     if ( fCurrentChannel == -1 ) {
 245         fCurrentChannel = iChannel;
 246         fCurrentSector  = icsector;
 247         fCurrentRow     = icRow;
 248     }
 249
 250     //process last pad if we change to a new one
 251     if ( iChannel != fCurrentChannel ){
 252         ProcessPad();
 253         fCurrentChannel = iChannel;
 254         fCurrentSector  = icsector;
 255         fCurrentRow     = icRow;
 256     }
 257
 258     //fill signals for current pad
 259     fPadSignal[icTimeBin]=csignal;
 260     if ( csignal > fMaxPadSignal ){
 261         fMaxPadSignal = csignal;
 262         fMaxTimeBin   = icTimeBin;
 263     }
 264     return 0;
 265 }
 266 //_____________________________________________________________________
 267 void AliTPCCalibSignal::ProcessPad() /*FOLD00*/
 268 {
 269     //
 270     //  Process data of current pad
 271     //
 272
 273     Float_t pedestal = 0;
 274
 275     if ( fBpedestal ){
 276         //!!!!!!! does not work like this
 277         //use pedestal database
 278         AliTPCCalROC *pedestalROC = 0x0;
 279
 280         //only load new pedestals if the sector has changed
 281         if ( fCurrentSector!=fLastSector ){
 282             pedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
 283             fLastSector=fCurrentSector;
 284         }
 285
 286         pedestal = pedestalROC->GetValue(fCurrentChannel);
 287
 288     } else {
 289
 290         //find pedestal for pad on the fly
 291         //using a few timebins before the signal
 292         Int_t pminus1 = 10, pminus2=5;
 293         Float_t sumN=0;
 294
 295         for (Int_t i=fMaxTimeBin-pminus1; i<fMaxTimeBin-pminus2+1; i++){
 296             if ( i>fFirstTimeBin && i<fLastTimeBin ){
 297                 pedestal+=fPadSignal[i];
 298                 sumN+=1.;
 299             }
 300         }
 301
 302         if ( sumN>0 ) pedestal/=sumN;
 303     }
 304
 305
 306
 307     //!!!! check borders
 308     //find signal mean and sigma
 309     Int_t tminus = 2, tplus=7;
 310     Double_t meanT=0, sigmaT=0, Qsum=0;
 311
 312
 313     for (Int_t i=fMaxTimeBin-tminus; i<fMaxTimeBin+tplus; i++){
 314         if ( i>=fFirstTimeBin && i<=fLastTimeBin ){
 315             Double_t val=fPadSignal[i]-pedestal;
 316             meanT+=val*(i+.5);      //+.5: center of the timebin
 317             sigmaT+=val*(i+.5)*(i+.5);
 318             Qsum+=val;
 319         }
 320     }
 321
 322
 323
 324     //!!!! What to do if Qsum == 0???
 325     //!!!! Should there be some threshold for max - pedestal and/or Qsum???
 326     //!!!! What if Qsum < 0
 327     //!!!! only fill time0 offset if Qsum > 0???
 328     if ( Qsum > 0 ){
 329         meanT/=Qsum;
 330         sigmaT/=Qsum;
 331         sigmaT = TMath::Sqrt(TMath::Abs(meanT*meanT - sigmaT));
 332
 333         //fill Time0 offset data for this event
 334         fVTime0Offset1[fCurrentSector]+=meanT;
 335         fVTime0Offset1Counter[fCurrentSector]++;
 336     } else {
 337         Qsum=0;
 338         meanT  = fLastTimeBinT0+1;               //put to overflow bin
 339         sigmaT = fLastTimeBinT0-fFirstTimeBinT0; //put to overflow bin
 340     }
 341
 342     //Fill Event T0 counter
 343     (*GetPadTimesEvent(fCurrentSector,kTRUE))[fCurrentChannel] = meanT;
 344
 345
 346     //Normalise Q to pad area of irocs
 347     Float_t norm = fParam->GetPadPitchWidth(0)*fParam->GetPadPitchLength(0,0)/(
 348         fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow));
 349
 350     //Fill Q histogram
 351     GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum*norm), fCurrentChannel );
 352
 353     //Fill RMS histogram
 354     GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
 355
 356
 357     //Fill debugging info
 358     if ( fDebugLevel>0 ){
 359         (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=pedestal;
 360         (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
 361         (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=Qsum;
 362     }
 363
 364     ResetPad();
 365 }
 366 //_____________________________________________________________________
 367 void AliTPCCalibSignal::EndEvent() /*FOLD00*/
 368 {
 369     //
 370     //  Process data of current pad
 371     //
 372     //check if last pad has allready been processed, if not do so
 373     if ( fMaxTimeBin>-1 ) ProcessPad();
 374
 375     //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
 376     for ( Int_t iSec = 0; iSec<72; iSec++ ){
 377         TVectorF *vTimes = GetPadTimesEvent(iSec);
 378         if ( !vTimes ) continue;
 379
 380         for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++ ){
 381             Float_t Time0 = fVTime0Offset1[iSec]/fVTime0Offset1Counter[iSec];
 382             Float_t Time  = (*vTimes)[iChannel];
 383
 384             GetHistoT0(iSec,kTRUE)->Fill( Time-Time0,iChannel );
 385
 386
 387             //Debug start
 388             if ( fDebugLevel>0 ){
 389                 if ( !fDebugStreamer ) {
 390                         //debug stream
 391                     TDirectory *backup = gDirectory;
 392                     fDebugStreamer = new TTreeSRedirector("deb2.root");
 393                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 394                 }
 395
 396                 Int_t row=0;
 397                 Int_t pad=0;
 398                 Int_t padc=0;
 399
 400                 Float_t Q   = (*GetPadQEvent(iSec))[iChannel];
 401                 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
 402
 403                 UInt_t channel=iChannel;
 404                 Int_t sector=iSec;
 405
 406                 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
 407                 pad = channel-fROC->GetRowIndexes(sector)[row];
 408                 padc = pad-(fROC->GetNPads(sector,row)/2);
 409
 410                 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
 411                                     Form("hSignalD%d.%d.%d",sector,row,pad),
 412                                     fLastTimeBin-fFirstTimeBin,
 413                                     fFirstTimeBin,fLastTimeBin);
 414                 h1->SetDirectory(0);
 415
 416                 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
 417                     h1->Fill(i,fPadSignal(i));
 418
 419                 (*fDebugStreamer) << "DataPad" <<
 420                     "Event=" << fEvent <<
 421                     "Sector="<< sector <<
 422                     "Row="   << row<<
 423                     "Pad="   << pad <<
 424                     "PadC="  << padc <<
 425                     "PadSec="<< channel <<
 426                     "Time0="  << Time0 <<
 427                     "Time="  << Time <<
 428                     "RMS="   << RMS <<
 429                     "Sum="   << Q <<
 430                     "hist.=" << h1 <<
 431                     "\n";
 432
 433                 delete h1;
 434             }
 435             //Debug end
 436
 437         }
 438     }
 439
 440 }
 441 //_____________________________________________________________________
 442 Bool_t AliTPCCalibSignal::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
 443 {
 444   //
 445   // Event Processing loop - AliTPCRawStream
 446   //
 447
 448   rawStream->SetOldRCUFormat(1);
 449
 450   ResetEvent();
 451
 452   Bool_t withInput = kFALSE;
 453
 454   while (rawStream->Next()) {
 455
 456       Int_t isector  = rawStream->GetSector();                       //  current sector
 457       Int_t iRow     = rawStream->GetRow();                          //  current row
 458       Int_t iPad     = rawStream->GetPad();                          //  current pad
 459       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 460       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 461
 462       Update(isector,iRow,iPad,iTimeBin,signal);
 463       withInput = kTRUE;
 464   }
 465
 466   if (withInput){
 467       EndEvent();
 468   }
 469
 470   return withInput;
 471 }
 472 //_____________________________________________________________________
 473 Bool_t AliTPCCalibSignal::ProcessEvent(AliRawReader *rawReader)
 474 {
 475   //
 476   //  Event processing loop - AliRawReader
 477   //
 478
 479
 480   AliTPCRawStream rawStream(rawReader);
 481
 482   rawReader->Select("TPC");
 483
 484   return ProcessEvent(&rawStream);
 485 }
 486 //_____________________________________________________________________
 487 Bool_t AliTPCCalibSignal::ProcessEvent(eventHeaderStruct *event)
 488 {
 489   //
 490   //  Event processing loop - date event
 491   //
 492     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 493     Bool_t result=ProcessEvent(rawReader);
 494     delete rawReader;
 495     return result;
 496
 497 }
 498 //_____________________________________________________________________
 499 TH2S* AliTPCCalibSignal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
 500                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 501                                   Char_t *type, Bool_t force)
 502 {
 503     //
 504     // return pointer to Q histogram
 505     // if force is true create a new histogram if it doesn't exist allready
 506     //
 507     if ( !force || arr->UncheckedAt(sector) )
 508         return (TH2S*)arr->UncheckedAt(sector);
 509
 510     // if we are forced and histogram doesn't yes exist create it
 511     Char_t name[255], title[255];
 512
 513     sprintf(name,"hCalib%s%.2d",type,sector);
 514     sprintf(title,"%s calibration histogram sector %.2d",type,sector);
 515
 516     // new histogram with Q calib information. One value for each pad!
 517     TH2S* hist = new TH2S(name,title,
 518                           nbinsY, ymin, ymax,
 519                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 520     hist->SetDirectory(0);
 521     arr->AddAt(hist,sector);
 522     return hist;
 523 }
 524 //_____________________________________________________________________
 525 TH2S* AliTPCCalibSignal::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
 526 {
 527     //
 528     // return pointer to T0 histogram
 529     // if force is true create a new histogram if it doesn't exist allready
 530     //
 531     TObjArray *arr = &fHistoT0Array;
 532     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 533 }
 534 //_____________________________________________________________________
 535 TH2S* AliTPCCalibSignal::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
 536 {
 537     //
 538     // return pointer to Q histogram
 539     // if force is true create a new histogram if it doesn't exist allready
 540     //
 541     TObjArray *arr = &fHistoQArray;
 542     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 543 }
 544 //_____________________________________________________________________
 545 TH2S* AliTPCCalibSignal::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
 546 {
 547     //
 548     // return pointer to Q histogram
 549     // if force is true create a new histogram if it doesn't exist allready
 550     //
 551     TObjArray *arr = &fHistoRMSArray;
 552     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 553 }
 554 //_____________________________________________________________________
 555 TVectorF* AliTPCCalibSignal::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
 556 {
 557     //
 558     // return pointer to Pad Info from 'arr' for the current event and sector
 559     // if force is true create it if it doesn't exist allready
 560     //
 561     if ( !force || arr->UncheckedAt(sector) )
 562         return (TVectorF*)arr->UncheckedAt(sector);
 563
 564     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 565     arr->AddAt(vect,sector);
 566     return vect;
 567 }
 568 //_____________________________________________________________________
 569 TVectorF* AliTPCCalibSignal::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
 570 {
 571     //
 572     // return pointer to Pad Times Array for the current event and sector
 573     // if force is true create it if it doesn't exist allready
 574     //
 575     TObjArray *arr = &fPadTimesArrayEvent;
 576     return GetPadInfoEvent(sector,arr,force);
 577 }
 578 //_____________________________________________________________________
 579 TVectorF* AliTPCCalibSignal::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
 580 {
 581     //
 582     // return pointer to Pad Q Array for the current event and sector
 583     // if force is true create it if it doesn't exist allready
 584     // for debugging purposes only
 585     //
 586
 587     TObjArray *arr = &fPadQArrayEvent;
 588     return GetPadInfoEvent(sector,arr,force);
 589 }
 590 //_____________________________________________________________________
 591 TVectorF* AliTPCCalibSignal::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
 592 {
 593     //
 594     // return pointer to Pad RMS Array for the current event and sector
 595     // if force is true create it if it doesn't exist allready
 596     // for debugging purposes only
 597     //
 598     TObjArray *arr = &fPadRMSArrayEvent;
 599     return GetPadInfoEvent(sector,arr,force);
 600 }
 601 //_____________________________________________________________________
 602 TVectorF* AliTPCCalibSignal::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
 603 {
 604     //
 605     // return pointer to Pad RMS Array for the current event and sector
 606     // if force is true create it if it doesn't exist allready
 607     // for debugging purposes only
 608     //
 609     TObjArray *arr = &fPadPedestalArrayEvent;
 610     return GetPadInfoEvent(sector,arr,force);
 611 }
 612 //_____________________________________________________________________
 613 AliTPCCalROC* AliTPCCalibSignal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
 614 {
 615     //
 616     // return pointer to ROC Calibration
 617     // if force is true create a new histogram if it doesn't exist allready
 618     //
 619     if ( !force || arr->UncheckedAt(sector) )
 620         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 621
 622     // if we are forced and histogram doesn't yes exist create it
 623
 624     // new AliTPCCalROC for T0 information. One value for each pad!
 625     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 626     //init values
 627     for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
 628         croc->SetValue(iChannel, 0);
 629     }
 630     arr->AddAt(croc,sector);
 631     return croc;
 632 }
 633 //_____________________________________________________________________
 634 AliTPCCalROC* AliTPCCalibSignal::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
 635 {
 636     //
 637     // return pointer to Carge ROC Calibration
 638     // if force is true create a new histogram if it doesn't exist allready
 639     //
 640     TObjArray *arr = &fCalRocArrayT0;
 641     return GetCalRoc(sector, arr, force);
 642 }
 643 //_____________________________________________________________________
 644 AliTPCCalROC* AliTPCCalibSignal::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
 645 {
 646     //
 647     // return pointer to T0 ROC Calibration
 648     // if force is true create a new histogram if it doesn't exist allready
 649     //
 650     TObjArray *arr = &fCalRocArrayQ;
 651     return GetCalRoc(sector, arr, force);
 652 }
 653 //_____________________________________________________________________
 654 AliTPCCalROC* AliTPCCalibSignal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
 655 {
 656     //
 657     // return pointer to signal width ROC Calibration
 658     // if force is true create a new histogram if it doesn't exist allready
 659     //
 660     TObjArray *arr = &fCalRocArrayRMS;
 661     return GetCalRoc(sector, arr, force);
 662 }
 663 //_____________________________________________________________________
 664 AliTPCCalROC* AliTPCCalibSignal::GetCalRocOutliers(Int_t sector, Bool_t force)
 665 {
 666     //
 667     // return pointer to Outliers
 668     // if force is true create a new histogram if it doesn't exist allready
 669     //
 670     TObjArray *arr = &fCalRocArrayOutliers;
 671     return GetCalRoc(sector, arr, force);
 672 }
 673 //_____________________________________________________________________
 674 void AliTPCCalibSignal::ResetEvent() /*FOLD00*/
 675 {
 676     //
 677     //  Reset global counters  -- Should be called before each event is processed
 678     //
 679     fLastSector=-1;
 680     fCurrentSector=-1;
 681     fCurrentRow=-1;
 682     fCurrentChannel=-1;
 683
 684     ResetPad();
 685
 686     fPadTimesArrayEvent.Delete();
 687     fPadQArrayEvent.Delete();
 688     fPadRMSArrayEvent.Delete();
 689     fPadPedestalArrayEvent.Delete();
 690
 691     for ( Int_t i=0; i<72; i++ ){
 692         fVTime0Offset1[i]=0;
 693         fVTime0Offset1Counter[i]=0;
 694     }
 695 }
 696 //_____________________________________________________________________
 697 void AliTPCCalibSignal::ResetPad() /*FOLD00*/
 698 {
 699     //
 700     //  Reset pad infos -- Should be called after a pad has been processed
 701     //
 702     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
 703         fPadSignal[i] = 0;
 704     fMaxTimeBin = -1;
 705     fMaxPadSignal = -1;
 706 }
 707 //_____________________________________________________________________
 708 void AliTPCCalibSignal::Analyse() /*FOLD00*/
 709 {
 710     //
 711     //  Calculate calibration constants
 712     //
 713
 714     TVectorD paramQ(3);
 715     TVectorD paramT0(3);
 716     TVectorD paramRMS(3);
 717     TMatrixD dummy(3,3);
 718
 719     for (Int_t iSec=0; iSec<72; iSec++){
 720         TH2S *hT0 = GetHistoT0(iSec);
 721         if (!hT0 ) continue;
 722
 723         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
 724         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
 725         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
 726         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
 727
 728         TH2S *hQ   = GetHistoQ(iSec);
 729         TH2S *hRMS = GetHistoRMS(iSec);
 730
 731         Short_t *array_hQ   = hQ->GetArray();
 732         Short_t *array_hT0  = hT0->GetArray();
 733         Short_t *array_hRMS = hRMS->GetArray();
 734
 735         UInt_t nChannels = fROC->GetNChannels(iSec);
 736
 737         //debug
 738         Int_t row=0;
 739         Int_t pad=0;
 740         Int_t padc=0;
 741         //! debug
 742
 743         for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
 744
 745
 746             Float_t cogTime0 = -1000;
 747             Float_t cogQ     = -1000;
 748             Float_t cogRMS   = -1000;
 749             Float_t cogOut   = 0;
 750
 751
 752             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
 753             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
 754             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
 755
 756 /*
 757             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
 758             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
 759             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
 760             cogQ     = paramQ[1];
 761             cogTime0 = paramT0[1];
 762             cogRMS   = paramRMS[1];
 763 */
 764             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
 765             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
 766             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
 767
 768
 769
 770             /*
 771             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
 772                 cogOut = 1;
 773                 cogTime0 = 0;
 774                 cogQ     = 0;
 775                 cogRMS   = 0;
 776             }
 777 */
 778             rocQ->SetValue(iChannel, cogQ*cogQ);
 779             rocT0->SetValue(iChannel, cogTime0);
 780             rocRMS->SetValue(iChannel, cogRMS);
 781             rocOut->SetValue(iChannel, cogOut);
 782
 783
 784             //debug
 785             if ( fDebugLevel > 0 ){
 786                 if ( !fDebugStreamer ) {
 787                         //debug stream
 788                     TDirectory *backup = gDirectory;
 789                     fDebugStreamer = new TTreeSRedirector("deb2.root");
 790                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 791                 }
 792
 793                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
 794                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
 795                 padc = pad-(fROC->GetNPads(iSec,row)/2);
 796
 797                 (*fDebugStreamer) << "DataEnd" <<
 798                     "Sector="  << iSec      <<
 799                     "Pad="     << pad       <<
 800                     "PadC="    << padc      <<
 801                     "Row="     << row       <<
 802                     "PadSec="  << iChannel   <<
 803                     "Q="       << cogQ      <<
 804                     "T0="      << cogTime0  <<
 805                     "RMS="     << cogRMS    <<
 806                     "\n";
 807             }
 808             //! debug
 809
 810         }
 811
 812     }
 813     delete fDebugStreamer;
 814     fDebugStreamer = 0x0;
 815 }
 816 //_____________________________________________________________________
 817 void AliTPCCalibSignal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
 818 {
 819     //
 820     //  Write class to file
 821     //
 822
 823     TString sDir(dir);
 824     TString option;
 825
 826     if ( append )
 827         option = "update";
 828     else
 829         option = "recreate";
 830
 831     TDirectory *backup = gDirectory;
 832     TFile f(filename,option.Data());
 833     f.cd();
 834     if ( !sDir.IsNull() ){
 835         f.mkdir(sDir.Data());
 836         f.cd(sDir);
 837     }
 838     this->Write();
 839     f.Close();
 840
 841     if ( backup ) backup->cd();
 842 }