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   ResetEvent();
 449
 450   Bool_t withInput = kFALSE;
 451
 452   while (rawStream->Next()) {
 453
 454       Int_t isector  = rawStream->GetSector();                       //  current sector
 455       Int_t iRow     = rawStream->GetRow();                          //  current row
 456       Int_t iPad     = rawStream->GetPad();                          //  current pad
 457       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 458       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 459
 460       Update(isector,iRow,iPad,iTimeBin,signal);
 461       withInput = kTRUE;
 462   }
 463
 464   if (withInput){
 465       EndEvent();
 466   }
 467
 468   return withInput;
 469 }
 470 //_____________________________________________________________________
 471 Bool_t AliTPCCalibSignal::ProcessEvent(AliRawReader *rawReader)
 472 {
 473   //
 474   //  Event processing loop - AliRawReader
 475   //
 476
 477
 478   AliTPCRawStream rawStream(rawReader);
 479
 480   rawReader->Select("TPC");
 481
 482   return ProcessEvent(&rawStream);
 483 }
 484 //_____________________________________________________________________
 485 Bool_t AliTPCCalibSignal::ProcessEvent(eventHeaderStruct *event)
 486 {
 487   //
 488   //  Event processing loop - date event
 489   //
 490     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 491     Bool_t result=ProcessEvent(rawReader);
 492     delete rawReader;
 493     return result;
 494
 495 }
 496 //_____________________________________________________________________
 497 TH2S* AliTPCCalibSignal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
 498                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 499                                   Char_t *type, Bool_t force)
 500 {
 501     //
 502     // return pointer to Q histogram
 503     // if force is true create a new histogram if it doesn't exist allready
 504     //
 505     if ( !force || arr->UncheckedAt(sector) )
 506         return (TH2S*)arr->UncheckedAt(sector);
 507
 508     // if we are forced and histogram doesn't yes exist create it
 509     Char_t name[255], title[255];
 510
 511     snprintf(name,255,"hCalib%s%.2d",type,sector);
 512     snprintf(title,255,"%s calibration histogram sector %.2d",type,sector);
 513
 514     // new histogram with Q calib information. One value for each pad!
 515     TH2S* hist = new TH2S(name,title,
 516                           nbinsY, ymin, ymax,
 517                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
 518     hist->SetDirectory(0);
 519     arr->AddAt(hist,sector);
 520     return hist;
 521 }
 522 //_____________________________________________________________________
 523 TH2S* AliTPCCalibSignal::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
 524 {
 525     //
 526     // return pointer to T0 histogram
 527     // if force is true create a new histogram if it doesn't exist allready
 528     //
 529     TObjArray *arr = &fHistoT0Array;
 530     return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
 531 }
 532 //_____________________________________________________________________
 533 TH2S* AliTPCCalibSignal::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
 534 {
 535     //
 536     // return pointer to Q histogram
 537     // if force is true create a new histogram if it doesn't exist allready
 538     //
 539     TObjArray *arr = &fHistoQArray;
 540     return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
 541 }
 542 //_____________________________________________________________________
 543 TH2S* AliTPCCalibSignal::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
 544 {
 545     //
 546     // return pointer to Q histogram
 547     // if force is true create a new histogram if it doesn't exist allready
 548     //
 549     TObjArray *arr = &fHistoRMSArray;
 550     return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
 551 }
 552 //_____________________________________________________________________
 553 TVectorF* AliTPCCalibSignal::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
 554 {
 555     //
 556     // return pointer to Pad Info from 'arr' for the current event and sector
 557     // if force is true create it if it doesn't exist allready
 558     //
 559     if ( !force || arr->UncheckedAt(sector) )
 560         return (TVectorF*)arr->UncheckedAt(sector);
 561
 562     TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
 563     arr->AddAt(vect,sector);
 564     return vect;
 565 }
 566 //_____________________________________________________________________
 567 TVectorF* AliTPCCalibSignal::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
 568 {
 569     //
 570     // return pointer to Pad Times Array for the current event and sector
 571     // if force is true create it if it doesn't exist allready
 572     //
 573     TObjArray *arr = &fPadTimesArrayEvent;
 574     return GetPadInfoEvent(sector,arr,force);
 575 }
 576 //_____________________________________________________________________
 577 TVectorF* AliTPCCalibSignal::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
 578 {
 579     //
 580     // return pointer to Pad Q Array for the current event and sector
 581     // if force is true create it if it doesn't exist allready
 582     // for debugging purposes only
 583     //
 584
 585     TObjArray *arr = &fPadQArrayEvent;
 586     return GetPadInfoEvent(sector,arr,force);
 587 }
 588 //_____________________________________________________________________
 589 TVectorF* AliTPCCalibSignal::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
 590 {
 591     //
 592     // return pointer to Pad RMS Array for the current event and sector
 593     // if force is true create it if it doesn't exist allready
 594     // for debugging purposes only
 595     //
 596     TObjArray *arr = &fPadRMSArrayEvent;
 597     return GetPadInfoEvent(sector,arr,force);
 598 }
 599 //_____________________________________________________________________
 600 TVectorF* AliTPCCalibSignal::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
 601 {
 602     //
 603     // return pointer to Pad RMS Array for the current event and sector
 604     // if force is true create it if it doesn't exist allready
 605     // for debugging purposes only
 606     //
 607     TObjArray *arr = &fPadPedestalArrayEvent;
 608     return GetPadInfoEvent(sector,arr,force);
 609 }
 610 //_____________________________________________________________________
 611 AliTPCCalROC* AliTPCCalibSignal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
 612 {
 613     //
 614     // return pointer to ROC Calibration
 615     // if force is true create a new histogram if it doesn't exist allready
 616     //
 617     if ( !force || arr->UncheckedAt(sector) )
 618         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 619
 620     // if we are forced and histogram doesn't yes exist create it
 621
 622     // new AliTPCCalROC for T0 information. One value for each pad!
 623     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 624     //init values
 625     for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
 626         croc->SetValue(iChannel, 0);
 627     }
 628     arr->AddAt(croc,sector);
 629     return croc;
 630 }
 631 //_____________________________________________________________________
 632 AliTPCCalROC* AliTPCCalibSignal::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
 633 {
 634     //
 635     // return pointer to Carge ROC Calibration
 636     // if force is true create a new histogram if it doesn't exist allready
 637     //
 638     TObjArray *arr = &fCalRocArrayT0;
 639     return GetCalRoc(sector, arr, force);
 640 }
 641 //_____________________________________________________________________
 642 AliTPCCalROC* AliTPCCalibSignal::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
 643 {
 644     //
 645     // return pointer to T0 ROC Calibration
 646     // if force is true create a new histogram if it doesn't exist allready
 647     //
 648     TObjArray *arr = &fCalRocArrayQ;
 649     return GetCalRoc(sector, arr, force);
 650 }
 651 //_____________________________________________________________________
 652 AliTPCCalROC* AliTPCCalibSignal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
 653 {
 654     //
 655     // return pointer to signal width ROC Calibration
 656     // if force is true create a new histogram if it doesn't exist allready
 657     //
 658     TObjArray *arr = &fCalRocArrayRMS;
 659     return GetCalRoc(sector, arr, force);
 660 }
 661 //_____________________________________________________________________
 662 AliTPCCalROC* AliTPCCalibSignal::GetCalRocOutliers(Int_t sector, Bool_t force)
 663 {
 664     //
 665     // return pointer to Outliers
 666     // if force is true create a new histogram if it doesn't exist allready
 667     //
 668     TObjArray *arr = &fCalRocArrayOutliers;
 669     return GetCalRoc(sector, arr, force);
 670 }
 671 //_____________________________________________________________________
 672 void AliTPCCalibSignal::ResetEvent() /*FOLD00*/
 673 {
 674     //
 675     //  Reset global counters  -- Should be called before each event is processed
 676     //
 677     fLastSector=-1;
 678     fCurrentSector=-1;
 679     fCurrentRow=-1;
 680     fCurrentChannel=-1;
 681
 682     ResetPad();
 683
 684     fPadTimesArrayEvent.Delete();
 685     fPadQArrayEvent.Delete();
 686     fPadRMSArrayEvent.Delete();
 687     fPadPedestalArrayEvent.Delete();
 688
 689     for ( Int_t i=0; i<72; i++ ){
 690         fVTime0Offset1[i]=0;
 691         fVTime0Offset1Counter[i]=0;
 692     }
 693 }
 694 //_____________________________________________________________________
 695 void AliTPCCalibSignal::ResetPad() /*FOLD00*/
 696 {
 697     //
 698     //  Reset pad infos -- Should be called after a pad has been processed
 699     //
 700     for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
 701         fPadSignal[i] = 0;
 702     fMaxTimeBin = -1;
 703     fMaxPadSignal = -1;
 704 }
 705 //_____________________________________________________________________
 706 void AliTPCCalibSignal::Analyse() /*FOLD00*/
 707 {
 708     //
 709     //  Calculate calibration constants
 710     //
 711
 712     TVectorD paramQ(3);
 713     TVectorD paramT0(3);
 714     TVectorD paramRMS(3);
 715     TMatrixD dummy(3,3);
 716
 717     for (Int_t iSec=0; iSec<72; iSec++){
 718         TH2S *hT0 = GetHistoT0(iSec);
 719         if (!hT0 ) continue;
 720
 721         AliTPCCalROC *rocQ   = GetCalRocQ  (iSec,kTRUE);
 722         AliTPCCalROC *rocT0  = GetCalRocT0 (iSec,kTRUE);
 723         AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
 724         AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
 725
 726         TH2S *hQ   = GetHistoQ(iSec);
 727         TH2S *hRMS = GetHistoRMS(iSec);
 728
 729         Short_t *array_hQ   = hQ->GetArray();
 730         Short_t *array_hT0  = hT0->GetArray();
 731         Short_t *array_hRMS = hRMS->GetArray();
 732
 733         UInt_t nChannels = fROC->GetNChannels(iSec);
 734
 735         //debug
 736         Int_t row=0;
 737         Int_t pad=0;
 738         Int_t padc=0;
 739         //! debug
 740
 741         for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
 742
 743
 744             Float_t cogTime0 = -1000;
 745             Float_t cogQ     = -1000;
 746             Float_t cogRMS   = -1000;
 747             Float_t cogOut   = 0;
 748
 749
 750             Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
 751             Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
 752             Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
 753
 754 /*
 755             AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
 756             AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
 757             AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
 758             cogQ     = paramQ[1];
 759             cogTime0 = paramT0[1];
 760             cogRMS   = paramRMS[1];
 761 */
 762             cogQ     = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
 763             cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
 764             cogRMS   = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
 765
 766
 767
 768             /*
 769             if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
 770                 cogOut = 1;
 771                 cogTime0 = 0;
 772                 cogQ     = 0;
 773                 cogRMS   = 0;
 774             }
 775 */
 776             rocQ->SetValue(iChannel, cogQ*cogQ);
 777             rocT0->SetValue(iChannel, cogTime0);
 778             rocRMS->SetValue(iChannel, cogRMS);
 779             rocOut->SetValue(iChannel, cogOut);
 780
 781
 782             //debug
 783             if ( fDebugLevel > 0 ){
 784                 if ( !fDebugStreamer ) {
 785                         //debug stream
 786                     TDirectory *backup = gDirectory;
 787                     fDebugStreamer = new TTreeSRedirector("deb2.root");
 788                     if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
 789                 }
 790
 791                 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
 792                 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
 793                 padc = pad-(fROC->GetNPads(iSec,row)/2);
 794
 795                 (*fDebugStreamer) << "DataEnd" <<
 796                     "Sector="  << iSec      <<
 797                     "Pad="     << pad       <<
 798                     "PadC="    << padc      <<
 799                     "Row="     << row       <<
 800                     "PadSec="  << iChannel   <<
 801                     "Q="       << cogQ      <<
 802                     "T0="      << cogTime0  <<
 803                     "RMS="     << cogRMS    <<
 804                     "\n";
 805             }
 806             //! debug
 807
 808         }
 809
 810     }
 811     delete fDebugStreamer;
 812     fDebugStreamer = 0x0;
 813 }
 814 //_____________________________________________________________________
 815 void AliTPCCalibSignal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
 816 {
 817     //
 818     //  Write class to file
 819     //
 820
 821     TString sDir(dir);
 822     TString option;
 823
 824     if ( append )
 825         option = "update";
 826     else
 827         option = "recreate";
 828
 829     TDirectory *backup = gDirectory;
 830     TFile f(filename,option.Data());
 831     f.cd();
 832     if ( !sDir.IsNull() ){
 833         f.mkdir(sDir.Data());
 834         f.cd(sDir);
 835     }
 836     this->Write();
 837     f.Close();
 838
 839     if ( backup ) backup->cd();
 840 }