TPC/AliTPCCalibPedestal.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 //          Implementation of the TPC pedestal and noise calibration
  19 //
  20 //   Origin: Jens Wiechula, Marian Ivanov   J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
  21 //
  22 //
  23 //-------------------------------------------------------
  24
  25
  26 /* $Id$ */
  27
  28 /*
  29  example: fill pedestal with gausschen noise
  30  AliTPCCalibPedestal ped;
  31  ped.TestEvent();
  32  ped.Analyse();
  33  //Draw output;
  34  TCanvas* c1 = new TCanvas;
  35  c1->Divide(1,2);
  36  c1->cd(1);
  37  ped.GetHistoPedestal(0)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
  38  ped.GetHistoPedestal(0)->Draw("colz");
  39  c1->cd(2);
  40  ped.GetHistoPedestal(36)->SetEntries(1); //needed in order for colz to work, reason: fast filling does not increase the entries counter
  41  ped.GetHistoPedestal(36)->Draw("colz");
  42  TCanvas* c2 = new TCanvas;
  43  c2->Divide(2,2);
  44  c2->cd(1);
  45  ped.GetCalRocPedestal(0)->Draw("colz");
  46  c2->cd(2);
  47  ped.GetCalRocRMS(0)->Draw("colz");
  48  c2->cd(3);
  49  ped.GetCalRocPedestal(36)->Draw("colz");
  50  c2->cd(4);
  51  ped.GetCalRocRMS(36)->Draw("colz");
  52
  53
  54 */
  55
  56 //Root includes
  57 #include <TObjArray.h>
  58 #include <TH1F.h>
  59 #include <TH2F.h>
  60 #include <TString.h>
  61 #include <TMath.h>
  62 #include <TF1.h>
  63 #include <TRandom.h>
  64 #include <TDirectory.h>
  65 #include <TFile.h>
  66 //AliRoot includes
  67 #include "AliRawReader.h"
  68 #include "AliRawReaderRoot.h"
  69 #include "AliRawReaderDate.h"
  70 #include "AliTPCRawStream.h"
  71 #include "AliTPCCalROC.h"
  72 #include "AliTPCROC.h"
  73 #include "AliMathBase.h"
  74 #include "TTreeStream.h"
  75
  76 //date
  77 #include "event.h"
  78
  79 //header file
  80 #include "AliTPCCalibPedestal.h"
  81
  82
  83 ClassImp(AliTPCCalibPedestal) /*FOLD00*/
  84
  85
  86
  87
  88
  89
  90
  91 AliTPCCalibPedestal::AliTPCCalibPedestal() : /*FOLD00*/
  92   TObject(),
  93   fFirstTimeBin(60),
  94   fLastTimeBin(1000),
  95   fAdcMin(1),
  96   fAdcMax(100),
  97   fROC(AliTPCROC::Instance()),
  98   fCalRocArrayPedestal(72),
  99   fCalRocArrayRMS(72),
 100   fHistoPedestalArray(72)
 101 {
 102     //
 103     // default constructor
 104     //
 105 }
 106 //_____________________________________________________________________
 107 AliTPCCalibPedestal::AliTPCCalibPedestal(const AliTPCCalibPedestal &ped) : /*FOLD00*/
 108   TObject(ped),
 109   fFirstTimeBin(ped.GetFirstTimeBin()),
 110   fLastTimeBin(ped.GetLastTimeBin()),
 111   fAdcMin(ped.GetAdcMin()),
 112   fAdcMax(ped.GetAdcMax()),
 113   fROC(AliTPCROC::Instance()),
 114   fCalRocArrayPedestal(72),
 115   fCalRocArrayRMS(72),
 116   fHistoPedestalArray(72)
 117 {
 118     //
 119     // copy constructor
 120     //
 121     for (Int_t iSec = 0; iSec < 72; iSec++){
 122         const AliTPCCalROC *calPed = (AliTPCCalROC*)ped.fCalRocArrayPedestal.UncheckedAt(iSec);
 123         const AliTPCCalROC *calRMS = (AliTPCCalROC*)ped.fCalRocArrayRMS.UncheckedAt(iSec);
 124         const TH2F         *hPed   = (TH2F*)ped.fHistoPedestalArray.UncheckedAt(iSec);
 125
 126         if ( calPed != 0x0 ) fCalRocArrayPedestal.AddAt(new AliTPCCalROC(*calPed), iSec);
 127         if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
 128
 129         if ( hPed != 0x0 ){
 130             TH2F *hNew = new TH2F(*hPed);
 131             hNew->SetDirectory(0);
 132             fHistoPedestalArray.AddAt(hNew,iSec);
 133         }
 134     }
 135 }
 136 //_____________________________________________________________________
 137 AliTPCCalibPedestal& AliTPCCalibPedestal::operator = (const  AliTPCCalibPedestal &source)
 138 {
 139   //
 140   // assignment operator
 141   //
 142   if (&source == this) return *this;
 143   new (this) AliTPCCalibPedestal(source);
 144
 145   return *this;
 146 }
 147 //_____________________________________________________________________
 148 AliTPCCalibPedestal::~AliTPCCalibPedestal() /*FOLD00*/
 149 {
 150   //
 151   // destructor
 152   //
 153   delete fROC;
 154 }
 155
 156
 157
 158
 159 //_____________________________________________________________________
 160 Int_t AliTPCCalibPedestal::Update(const Int_t icsector, /*FOLD00*/
 161                                 const Int_t icRow,
 162                                 const Int_t icPad,
 163                                 const Int_t icTimeBin,
 164                                 const Float_t csignal)
 165 {
 166     //
 167     // Signal filling methode
 168     //
 169     if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin)   ) return 0;
 170
 171     Int_t iChannel  = fROC->GetRowIndexes(icsector)[icRow]+icPad; //  global pad position in sector
 172
 173     // fast filling methode.
 174     // Attention: the entry counter of the histogram is not increased
 175     //            this means that e.g. the colz draw option gives an empty plot
 176     Int_t bin = 0;
 177     if ( !(((Int_t)csignal>fAdcMax ) || ((Int_t)csignal<fAdcMin)) )
 178         bin = (iChannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
 179
 180     GetHistoPedestal(icsector,kTRUE)->GetArray()[bin]++;
 181
 182     return 0;
 183 }
 184 //_____________________________________________________________________
 185 Bool_t AliTPCCalibPedestal::ProcessEvent(AliTPCRawStream *rawStream)
 186 {
 187   //
 188   // Event Processing loop - AliTPCRawStream
 189   //
 190
 191   rawStream->SetOldRCUFormat(1);
 192
 193   Bool_t withInput = kFALSE;
 194
 195   while (rawStream->Next()) {
 196
 197       Int_t isector  = rawStream->GetSector();                       //  current sector
 198       Int_t iRow     = rawStream->GetRow();                          //  current row
 199       Int_t iPad     = rawStream->GetPad();                          //  current pad
 200       Int_t iTimeBin = rawStream->GetTime();                         //  current time bin
 201       Float_t signal = rawStream->GetSignal();                       //  current ADC signal
 202
 203       Update(isector,iRow,iPad,iTimeBin,signal);
 204       withInput = kTRUE;
 205   }
 206
 207   return withInput;
 208 }
 209 //_____________________________________________________________________
 210 Bool_t AliTPCCalibPedestal::ProcessEvent(AliRawReader *rawReader)
 211 {
 212   //
 213   //  Event processing loop - AliRawReader
 214   //
 215
 216
 217   AliTPCRawStream rawStream(rawReader);
 218
 219   rawReader->Select("TPC");
 220
 221   return ProcessEvent(&rawStream);
 222 }
 223 //_____________________________________________________________________
 224 Bool_t AliTPCCalibPedestal::ProcessEvent(eventHeaderStruct *event)
 225 {
 226   //
 227   //  process date event
 228   //
 229     AliRawReader *rawReader = new AliRawReaderDate((void*)event);
 230     Bool_t result=ProcessEvent(rawReader);
 231     delete rawReader;
 232     return result;
 233 }
 234 //_____________________________________________________________________
 235 Bool_t AliTPCCalibPedestal::TestEvent() /*FOLD00*/
 236 {
 237   //
 238   //  Test event loop
 239   // fill one oroc and one iroc with random gaus
 240   //
 241
 242     gRandom->SetSeed(0);
 243
 244     for (UInt_t iSec=0; iSec<72; iSec++){
 245         if (iSec%36>0) continue;
 246         for (UInt_t iRow=0; iRow < fROC->GetNRows(iSec); iRow++){
 247             for (UInt_t iPad=0; iPad < fROC->GetNPads(iSec,iRow); iPad++){
 248                 for (UInt_t iTimeBin=0; iTimeBin<1024; iTimeBin++){
 249                     Float_t signal=(Int_t)(iRow+3+gRandom->Gaus(0,.7));
 250                     if ( signal>0 )Update(iSec,iRow,iPad,iTimeBin,signal);
 251                 }
 252             }
 253         }
 254     }
 255     return kTRUE;
 256 }
 257 //_____________________________________________________________________
 258 TH2F* AliTPCCalibPedestal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
 259                                   Int_t nbinsY, Float_t ymin, Float_t ymax,
 260                                   Char_t *type, Bool_t force)
 261 {
 262     //
 263     // return pointer to Q histogram
 264     // if force is true create a new histogram if it doesn't exist allready
 265     //
 266     if ( !force || arr->UncheckedAt(sector) )
 267         return (TH2F*)arr->UncheckedAt(sector);
 268
 269     // if we are forced and histogram doesn't yes exist create it
 270     Char_t name[255], title[255];
 271
 272     sprintf(name,"hCalib%s%.2d",type,sector);
 273     sprintf(title,"%s calibration histogram sector %.2d;ADC channel;Channel (pad)",type,sector);
 274
 275     // new histogram with Q calib information. One value for each pad!
 276     TH2F* hist = new TH2F(name,title,
 277                           nbinsY, ymin, ymax,
 278                           fROC->GetNChannels(sector),0,fROC->GetNChannels(sector)
 279                          );
 280     hist->SetDirectory(0);
 281     arr->AddAt(hist,sector);
 282     return hist;
 283 }
 284 //_____________________________________________________________________
 285 TH2F* AliTPCCalibPedestal::GetHistoPedestal(Int_t sector, Bool_t force) /*FOLD00*/
 286 {
 287     //
 288     // return pointer to T0 histogram
 289     // if force is true create a new histogram if it doesn't exist allready
 290     //
 291     TObjArray *arr = &fHistoPedestalArray;
 292     return GetHisto(sector, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
 293 }
 294 //_____________________________________________________________________
 295 AliTPCCalROC* AliTPCCalibPedestal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
 296 {
 297     //
 298     // return pointer to ROC Calibration
 299     // if force is true create a new histogram if it doesn't exist allready
 300     //
 301     if ( !force || arr->UncheckedAt(sector) )
 302         return (AliTPCCalROC*)arr->UncheckedAt(sector);
 303
 304     // if we are forced and histogram doesn't yes exist create it
 305
 306     // new AliTPCCalROC for T0 information. One value for each pad!
 307     AliTPCCalROC *croc = new AliTPCCalROC(sector);
 308     //init values
 309     for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
 310         croc->SetValue(iChannel, 0);
 311     }
 312     arr->AddAt(croc,sector);
 313     return croc;
 314 }
 315 //_____________________________________________________________________
 316 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocPedestal(Int_t sector, Bool_t force) /*FOLD00*/
 317 {
 318     //
 319     // return pointer to Carge ROC Calibration
 320     // if force is true create a new histogram if it doesn't exist allready
 321     //
 322     TObjArray *arr = &fCalRocArrayPedestal;
 323     return GetCalRoc(sector, arr, force);
 324 }
 325 //_____________________________________________________________________
 326 AliTPCCalROC* AliTPCCalibPedestal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
 327 {
 328     //
 329     // return pointer to signal width ROC Calibration
 330     // if force is true create a new histogram if it doesn't exist allready
 331     //
 332     TObjArray *arr = &fCalRocArrayRMS;
 333     return GetCalRoc(sector, arr, force);
 334 }
 335 //_____________________________________________________________________
 336 void AliTPCCalibPedestal::Analyse() /*FOLD00*/
 337 {
 338     //
 339     //  Calculate calibration constants
 340     //
 341
 342     Int_t nbinsAdc = fAdcMax-fAdcMin;
 343
 344     TVectorD param(3);
 345     TMatrixD dummy(3,3);
 346
 347     Float_t *array_hP=0;
 348
 349
 350     for (Int_t iSec=0; iSec<72; iSec++){
 351         TH2F *hP = GetHistoPedestal(iSec);
 352         if ( !hP ) continue;
 353
 354         AliTPCCalROC *rocPedestal = GetCalRocPedestal(iSec,kTRUE);
 355         AliTPCCalROC *rocRMS      = GetCalRocRMS(iSec,kTRUE);
 356
 357         array_hP = hP->GetArray();
 358         UInt_t nChannels = fROC->GetNChannels(iSec);
 359
 360         for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
 361             Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
 362             Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,&param,&dummy);
 363             // if the fitting failed set noise and pedestal to 0
 364             if ( ret == -4 ) {
 365                 param[1]=0;
 366                 param[2]=0;
 367             }
 368             rocPedestal->SetValue(iChannel,param[1]);
 369             rocRMS->SetValue(iChannel,param[2]);
 370         }
 371     }
 372 }
 373 //_____________________________________________________________________
 374 void AliTPCCalibPedestal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
 375 {
 376     //
 377     //  Write class to file
 378     //
 379
 380     TString sDir(dir);
 381     TString option;
 382
 383     if ( append )
 384         option = "update";
 385     else
 386         option = "recreate";
 387
 388     TDirectory *backup = gDirectory;
 389     TFile f(filename,option.Data());
 390     f.cd();
 391     if ( !sDir.IsNull() ){
 392         f.mkdir(sDir.Data());
 393         f.cd(sDir);
 394     }
 395     this->Write();
 396     f.Close();
 397
 398     if ( backup ) backup->cd();
 399 }