FMD/AliFMDPedestalDA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2004, 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 /** @file    AliFMDPedestalDA.cxx
  17     @author  Hans Hjersing Dalsgaard <canute@nbi.dk>
  18     @date    Mon Mar 10 09:46:05 2008
  19     @brief   Derived class for the pedestal detector algorithm.
  20 */
  21 //
  22 // This class implements the virtual functions of the AliFMDBaseDA
  23 // class.  The most important of these functions, FillChannels(..) and
  24 // Analyse(...) collect and analyse the data of each channel. The
  25 // resulting pedestal and noise values are written to a comma
  26 // separated values (csv) file on the go. The csv files produced in
  27 // this way are the basic input to the AliFMDPreprocessor.
  28 //
  29
  30 #include "AliFMDPedestalDA.h"
  31 #include "AliFMDAltroMapping.h"
  32 #include "AliFMDParameters.h"
  33 #include "AliFMDCalibPedestal.h"
  34 #include "AliFMDDigit.h"
  35 #include "AliLog.h"
  36 #include <iostream>
  37 #include <fstream>
  38 #include <iomanip>
  39 #include <TFile.h>
  40 #include <TF1.h>
  41 #include <TObject.h>
  42 #include <TMath.h>
  43 #include <TSystem.h>
  44 #include <TDatime.h>
  45 #include <TH2.h>
  46
  47 //_____________________________________________________________________
  48 ClassImp(AliFMDPedestalDA)
  49
  50 //_____________________________________________________________________
  51 AliFMDPedestalDA::AliFMDPedestalDA()
  52   : AliFMDBaseDA(),
  53     fCurrentChannel(1),
  54     fPedSummary("PedestalSummary","pedestals",51200,0,51200),
  55     fNoiseSummary("NoiseSummary","noise",51200,0,51200),
  56     fZSfileFMD1(),
  57     fZSfileFMD2(),
  58     fZSfileFMD3(),
  59     fMinTimebin(3 * 4 * 3 * 16), // 3 ddls, 4 FECs, 3 Altros, 16 channels
  60     fMaxTimebin(3 * 4 * 3 * 16), // 3 ddls, 4 FECs, 3 Altros, 16 channels
  61     fSummaryFMD1i(0),
  62     fSummaryFMD2i(0),
  63     fSummaryFMD2o(0),
  64     fSummaryFMD3i(0),
  65     fSummaryFMD3o(0)
  66 {
  67   // Default constructor
  68   Rotate("peds.csv", 3);
  69   fOutputFile.open("peds.csv");
  70   Rotate("ddl3072.csv", 10);
  71   fZSfileFMD1.open("ddl3072.csv");
  72   Rotate("ddl3073.csv", 10);
  73   fZSfileFMD2.open("ddl3073.csv");
  74   Rotate("ddl3074.csv", 10);
  75   fZSfileFMD3.open("ddl3074.csv");
  76 }
  77
  78 //_____________________________________________________________________
  79 AliFMDPedestalDA::AliFMDPedestalDA(const AliFMDPedestalDA & pedDA) :
  80   AliFMDBaseDA(pedDA),
  81   fCurrentChannel(1),
  82   fPedSummary("PedestalSummary","pedestals",51200,0,51200),
  83   fNoiseSummary("NoiseSummary","noise",51200,0,51200),
  84   fZSfileFMD1(),
  85   fZSfileFMD2(),
  86   fZSfileFMD3(),
  87   fMinTimebin(pedDA.fMinTimebin),
  88   fMaxTimebin(pedDA.fMaxTimebin),
  89   fSummaryFMD1i(pedDA.fSummaryFMD1i),
  90   fSummaryFMD2i(pedDA.fSummaryFMD2i),
  91   fSummaryFMD2o(pedDA.fSummaryFMD2o),
  92   fSummaryFMD3i(pedDA.fSummaryFMD3i),
  93   fSummaryFMD3o(pedDA.fSummaryFMD3o)
  94 {
  95   // Copy constructor
  96 }
  97
  98 //_____________________________________________________________________
  99 AliFMDPedestalDA::~AliFMDPedestalDA()
 100 {
 101   // Destructor.
 102 }
 103
 104 //_____________________________________________________________________
 105 void AliFMDPedestalDA::Init()
 106 {
 107   // Initialise
 108   SetRequiredEvents(1000);
 109   fMinTimebin.Reset(1024);
 110   fMaxTimebin.Reset(-1);
 111
 112
 113 }
 114
 115 //_____________________________________________________________________
 116 void AliFMDPedestalDA::AddChannelContainer(TObjArray* sectorArray,
 117                                            UShort_t det,
 118                                            Char_t   ring,
 119                                            UShort_t sec,
 120                                            UShort_t strip)
 121 {
 122   // Add a channel to the containers.
 123   //
 124   // Parameters:
 125   //     sectorArray  Array of sectors
 126   //     det          Detector
 127   //     ring         Ring
 128   //     sec          Sector
 129   //     strip        Strip
 130   AliFMDParameters* pars        = AliFMDParameters::Instance();
 131   UInt_t            samples     = pars->GetSampleRate(det, ring, sec, strip);
 132   TObjArray*        sampleArray = new TObjArray(samples);
 133   sampleArray->SetOwner();
 134   for (UInt_t sample = 0; sample < samples; sample++) {
 135     TString name(Form("FMD%d%c[%02d,%03d]_%d", det,ring,sec,strip,sample));
 136     TH1S* hSample = new TH1S(name.Data(),name.Data(), 1024,-.5,1023.5);
 137     hSample->SetXTitle("ADC");
 138     hSample->SetYTitle("Events");
 139     hSample->SetDirectory(0);
 140     sampleArray->AddAt(hSample, sample);
 141   }
 142   sectorArray->AddAtAndExpand(sampleArray, strip);
 143 }
 144
 145 //_____________________________________________________________________
 146 void AliFMDPedestalDA::FillChannels(AliFMDDigit* digit)
 147 {
 148   // Fill ADC values from a digit into the corresponding histogram.
 149   //
 150   // Parameters:
 151   //    digit Digit to fill ADC values for.
 152   UShort_t det   = digit->Detector();
 153   Char_t   ring  = digit->Ring();
 154   UShort_t sec   = digit->Sector();
 155   UShort_t strip = digit->Strip();
 156
 157   AliFMDParameters* pars     = AliFMDParameters::Instance();
 158   UInt_t            samples  = pars->GetSampleRate(det, ring, sec, strip);
 159   for (UInt_t sample = 0; sample < samples; sample++) {
 160     TH1S* hSample = GetChannel(det, ring, sec, strip, sample);
 161     hSample->Fill(digit->Count(sample));
 162   }
 163
 164 }
 165
 166 //_____________________________________________________________________
 167 void AliFMDPedestalDA::MakeSummary(UShort_t det, Char_t ring)
 168 {
 169   //Create summary hists for FMD pedestals
 170   std::cout << "Making summary for FMD" << det << ring << " ..." << std::endl;
 171   switch (det) {
 172   case 1:
 173     fSummaryFMD1i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 174     break;
 175   case 2:
 176     switch (ring) {
 177     case 'I': case 'i':
 178       fSummaryFMD2i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 179       break;
 180     case 'O': case 'o':
 181       fSummaryFMD2o = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 182       break;
 183     }
 184     break;
 185   case 3:
 186     switch (ring) {
 187     case 'I': case 'i':
 188       fSummaryFMD3i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 189       break;
 190     case 'O': case 'o':
 191       fSummaryFMD3o = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 192       break;
 193     }
 194     break;
 195   }
 196 }
 197
 198 //_____________________________________________________________________
 199 void AliFMDPedestalDA::Analyse(UShort_t det,
 200                                Char_t   ring,
 201                                UShort_t sec,
 202                                UShort_t strip)
 203 {
 204   // Analyse a strip.  That is, compute the mean and spread of the ADC
 205   // spectra for all strips.  Also output on files the values.
 206   //
 207   // Parameters:
 208   //     det   Detector
 209   //     ring  Ring
 210   //     sec   Sector
 211   //     strip Strip.
 212   AliFMDParameters* pars     = AliFMDParameters::Instance();
 213   TH2* summary = 0;
 214   switch (det) {
 215   case 1: summary = fSummaryFMD1i; break;
 216   case 2:
 217     switch (ring) {
 218     case 'I':  summary = fSummaryFMD2i; break;
 219     case 'O':  summary = fSummaryFMD2o; break;
 220     }
 221     break;
 222   case 3:
 223     switch (ring) {
 224     case 'I':  summary = fSummaryFMD3i; break;
 225     case 'O':  summary = fSummaryFMD3o; break;
 226     }
 227     break;
 228   }
 229   static bool first = true;
 230   if (summary && first) {
 231     std::cout << "Filling summary " << summary->GetName() << std::endl;
 232     first = false;
 233   }
 234
 235   // Float_t           factor   = pars->GetPedestalFactor();
 236   UInt_t            samples  = pars->GetSampleRate(det, ring, sec, strip);
 237   for (UShort_t sample = 0; sample < samples; sample++) {
 238
 239     TH1S* hChannel = GetChannel(det, ring, sec, strip,sample);
 240     if(hChannel->GetEntries() == 0) {
 241       //AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
 242       //            det,ring,sec,strip));
 243       return;
 244     }
 245
 246     AliDebug(50, Form("Fitting FMD%d%c[%02d,%03d] with %d entries",
 247                       det,ring,sec,strip, int(hChannel->GetEntries())));
 248     TF1 fitFunc("fitFunc","gausn",0,300);
 249     fitFunc.SetParameters(100,100,1);
 250     hChannel->Fit("fitFunc","Q0+","",10,200);
 251
 252     Float_t mean = hChannel->GetMean();
 253     Float_t rms  = hChannel->GetRMS();
 254
 255
 256
 257     hChannel->GetXaxis()->SetRangeUser(mean-5*rms,mean+5*rms);
 258
 259     mean = hChannel->GetMean();
 260     rms  = hChannel->GetRMS();
 261
 262
 263     UShort_t ddl, board, altro, channel;
 264     UShort_t timebin;
 265
 266     pars->Detector2Hardware(det,ring,sec,strip,sample,
 267                             ddl,board,altro,channel,timebin);
 268     Int_t idx = HWIndex(ddl, board, altro, channel);
 269     if (idx >= 0) {
 270       fMinTimebin[idx] = TMath::Min(Short_t(timebin),   fMinTimebin[idx]);
 271       fMaxTimebin[idx] = TMath::Max(Short_t(timebin+1), fMaxTimebin[idx]);
 272     }
 273
 274     std::ostream* zsFile = 0;
 275     switch(det) {
 276     case 1:  zsFile = &fZSfileFMD1; break;
 277     case 2:  zsFile = &fZSfileFMD2; break;
 278     case 3:  zsFile = &fZSfileFMD3; break;
 279     default: AliWarning("Unknown sample!"); break;
 280
 281     }
 282     *zsFile << board   << ','
 283             << altro   << ','
 284             << channel << ','
 285             << timebin << ','
 286             << mean    << ','
 287             << rms     << "\n";
 288
 289     Float_t chi2ndf = 0;
 290
 291
 292     if(fitFunc.GetNDF())
 293       chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
 294
 295
 296     Int_t sampleToWrite = 2;
 297     if      (samples == 2) sampleToWrite = 1;
 298     else if (samples <  2) sampleToWrite = 0;
 299
 300     if(sample != sampleToWrite) continue;
 301
 302
 303     fOutputFile << det                         << ','
 304                 << ring                        << ','
 305                 << sec                         << ','
 306                 << strip                       << ','
 307                 << mean                        << ','
 308                 << rms                         << ','
 309                 << fitFunc.GetParameter(1)     << ','
 310                 << fitFunc.GetParameter(2)     << ','
 311                 << chi2ndf                     <<"\n";
 312
 313     if (summary) {
 314       Int_t bin = summary->FindBin(sec, strip);
 315       summary->SetBinContent(bin, mean);
 316       summary->SetBinError(bin, rms);
 317     }
 318
 319     if(fSaveHistograms  ) {
 320       gDirectory->cd(GetSectorPath(det, ring, sec, kTRUE));
 321       TH1F* sumPed   = dynamic_cast<TH1F*>(gDirectory->Get("Pedestals"));
 322       TH1F* sumNoise = dynamic_cast<TH1F*>(gDirectory->Get("Noise"));
 323       Int_t nStr = (ring == 'I' ? 512 : 256);
 324       if (!sumPed) {
 325         sumPed = new TH1F("Pedestals",
 326                           Form("Summary of pedestals in FMD%d%c[%02d]",
 327                                det, ring, sec),
 328                           nStr, -.5, nStr-.5);
 329         sumPed->SetXTitle("Strip");
 330         sumPed->SetYTitle("Pedestal [ADC]");
 331         sumPed->SetDirectory(gDirectory);
 332       }
 333       if (!sumNoise) {
 334         sumNoise = new TH1F("Noise",
 335                             Form("Summary of noise in FMD%d%c[%02d]",
 336                                  det, ring, sec),
 337                             nStr, -.5, nStr-.5);
 338         sumNoise->SetXTitle("Strip");
 339         sumNoise->SetYTitle("Noise [ADC]");
 340
 341         sumNoise->SetDirectory(gDirectory);
 342       }
 343       sumPed->SetBinContent(strip+1, mean);
 344       sumPed->SetBinError(strip+1, rms);
 345       sumNoise->SetBinContent(strip+1, rms);
 346
 347       if(sumNoise->GetEntries() == nStr)
 348         sumNoise->Write(sumNoise->GetName(),TObject::kOverwrite);
 349       if(sumPed->GetEntries() == nStr)
 350         sumPed->Write(sumPed->GetName(),TObject::kOverwrite);
 351
 352       fPedSummary.SetBinContent(fCurrentChannel,mean);
 353
 354       fNoiseSummary.SetBinContent(fCurrentChannel,rms);
 355       fCurrentChannel++;
 356
 357       gDirectory->cd(GetStripPath(det, ring, sec, strip, kTRUE));
 358       hChannel->GetXaxis()->SetRange(1,1024);
 359
 360       hChannel->Write();
 361     }
 362   }
 363 }
 364
 365 //_____________________________________________________________________
 366 void AliFMDPedestalDA::Terminate(TFile* diagFile)
 367 {
 368   // Called at the end of a job.  Fills in missing time-bins and
 369   // closes output files
 370   if(fSaveHistograms) {
 371     diagFile->cd();
 372
 373     fPedSummary.Write();
 374     fNoiseSummary.Write();
 375   }
 376   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
 377   for (Int_t i = 0; i < 3; i++) {
 378     std::ofstream& out = (i == 0 ? fZSfileFMD1 :
 379                           i == 1 ? fZSfileFMD2 :
 380                           fZSfileFMD3);
 381     if (out.is_open() && fSeenDetectors[i]) {
 382       FillinTimebins(out, map->Detector2DDL(i+1));
 383     }
 384     if (!fSeenDetectors[i]) {
 385       TString n(Form("ddl%d.csv",3072+map->Detector2DDL(i+1)));
 386       gSystem->Unlink(n.Data());
 387     }
 388   }
 389
 390 }
 391
 392 //_____________________________________________________________________
 393 void AliFMDPedestalDA::FillinTimebins(std::ofstream& out, UShort_t /*ddl*/)
 394 {
 395   //
 396   // Fill missing timebins
 397   //
 398 #if 0
 399   unsigned short  boards[] = { 0x0, 0x1, 0x10, 0x11, 0xFFFF };
 400   unsigned short* board    = boards;
 401   while ((*boards) != 0xFFFF) {
 402     for (UShort_t altro = 0; altro < 3; altro++) {
 403       for (UShort_t channel = 0; channel < 16; channel++) {
 404         Int_t idx = HWIndex(ddl, *board, altro, channel);
 405         if (idx < 0) {
 406           AliWarning(Form("Invalid index for %4d/0x%02x/0x%x/0x%x: %d",
 407                           ddl, *board, altro, channel, idx));
 408           continue;
 409         }
 410         Short_t min = fMinTimebin[idx];
 411         Short_t max = fMaxTimebin[idx];
 412
 413         // Channel not seen at all.
 414         if (min > 1023 || max < 0) continue;
 415
 416         out << "# Extra timebins for 0x" << std::hex
 417             << board << ',' << altro << ',' << channel
 418             << " got time-bins " << min << " to " << max-1
 419             << std::dec << std::endl;
 420
 421         for (UShort_t t = 15; t < min; t++)
 422           // Write a phony line
 423           out << board << "," << altro << "," << channel << ","
 424               << t << "," << 1023 << "," << 0 << std::endl;
 425
 426         for (UShort_t t = max; t < 1024; t++)
 427           // Write a phony line
 428           out << board << "," << altro << "," << channel << ","
 429               << t << "," << 1023 << "," << 0 << std::endl;
 430       } // channel loop
 431     } // altro loop
 432   } // board loop
 433   // Write trailer, and close
 434 #endif
 435   out.write("# EOF\n", 6);
 436   out.close();
 437 }
 438
 439 //_____________________________________________________________________
 440 void AliFMDPedestalDA::WriteHeaderToFile()
 441 {
 442   //
 443   // Write headers to output files
 444   //
 445   AliFMDParameters* pars       = AliFMDParameters::Instance();
 446   fOutputFile.write(Form("# %s \n",pars->GetPedestalShuttleID()),13);
 447   TDatime now;
 448   fOutputFile << "# This file created from run # " << fRunno
 449               << " @ " << now.AsString() << std::endl;
 450   fOutputFile.write("# Detector, "
 451                     "Ring, "
 452                     "Sector, "
 453                     "Strip, "
 454                     "Pedestal, "
 455                     "Noise, "
 456                     "Mu, "
 457                     "Sigma, "
 458                     "Chi2/NDF \n", 71);
 459
 460   std::ostream* zss[] = { &fZSfileFMD1, &fZSfileFMD2, &fZSfileFMD3, 0 };
 461   for (size_t i = 0; i < 3; i++)  {
 462     *(zss[i]) << "# FMD " << (i+1) << " pedestals \n"
 463               << "# board, "
 464               << "altro, "
 465               << "channel, "
 466               << "timebin, "
 467               << "pedestal, "
 468               << "noise\n";
 469     *(zss[i]) << "# This file created from run # " << fRunno
 470               << " @ " << now.AsString() << std::endl;
 471   }
 472 }
 473
 474 //_____________________________________________________________________
 475 TH1S* AliFMDPedestalDA::GetChannel(UShort_t det,
 476                                    Char_t   ring,
 477                                    UShort_t sec,
 478                                    UShort_t strip,
 479                                    UInt_t   sample)
 480 {
 481   // Get the histogram corresponding to a strip sample.
 482   //
 483   // Parameters:
 484   //     det    Detector
 485   //     ring   Ring
 486   //     sec    Sector
 487   //     strip  Strip
 488   //     sample Sample
 489   //
 490   // Return:
 491   //     ADC spectra of a strip.
 492   UShort_t   iring       = (ring == 'O' ? 0 : 1);
 493   TObjArray* detArray    = static_cast<TObjArray*>(fDetectorArray.At(det));
 494   TObjArray* ringArray   = static_cast<TObjArray*>(detArray->At(iring));
 495   TObjArray* secArray    = static_cast<TObjArray*>(ringArray->At(sec));
 496   TObjArray* sampleArray = static_cast<TObjArray*>(secArray->At(strip));
 497   TH1S*      hSample     = static_cast<TH1S*>(sampleArray->At(sample));
 498   return hSample;
 499
 500 }
 501
 502 //_____________________________________________________________________
 503 //
 504 //EOF
 505 //