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   fDiagnosticsFilename = "diagnosticsPedestal.root";
  77 }
  78
  79 //_____________________________________________________________________
  80 AliFMDPedestalDA::AliFMDPedestalDA(const AliFMDPedestalDA & pedDA) :
  81   AliFMDBaseDA(pedDA),
  82   fCurrentChannel(1),
  83   fPedSummary("PedestalSummary","pedestals",51200,0,51200),
  84   fNoiseSummary("NoiseSummary","noise",51200,0,51200),
  85   fZSfileFMD1(),
  86   fZSfileFMD2(),
  87   fZSfileFMD3(),
  88   fMinTimebin(pedDA.fMinTimebin),
  89   fMaxTimebin(pedDA.fMaxTimebin),
  90   fSummaryFMD1i(pedDA.fSummaryFMD1i),
  91   fSummaryFMD2i(pedDA.fSummaryFMD2i),
  92   fSummaryFMD2o(pedDA.fSummaryFMD2o),
  93   fSummaryFMD3i(pedDA.fSummaryFMD3i),
  94   fSummaryFMD3o(pedDA.fSummaryFMD3o)
  95 {
  96   // Copy constructor
  97 }
  98
  99 //_____________________________________________________________________
 100 AliFMDPedestalDA::~AliFMDPedestalDA()
 101 {
 102   // Destructor.
 103 }
 104
 105 //_____________________________________________________________________
 106 void AliFMDPedestalDA::Init()
 107 {
 108   // Initialise
 109   SetRequiredEvents(1000);
 110   fMinTimebin.Reset(1024);
 111   fMaxTimebin.Reset(-1);
 112
 113
 114 }
 115
 116 //_____________________________________________________________________
 117 void AliFMDPedestalDA::AddChannelContainer(TObjArray* sampleArray,
 118                                            UShort_t det,
 119                                            Char_t   ring,
 120                                            UShort_t sec,
 121                                            UShort_t strip)
 122 {
 123   // Add a channel to the containers.
 124   //
 125   // Parameters:
 126   //     sectorArray  Array of sectors
 127   //     det          Detector
 128   //     ring         Ring
 129   //     sec          Sector
 130   //     strip        Strip
 131   AliFMDParameters* pars        = AliFMDParameters::Instance();
 132   UInt_t            samples     = pars->GetSampleRate(det, ring, sec, strip);
 133   for (UInt_t sample = 0; sample < samples; sample++) {
 134     TString name(Form("FMD%d%c[%02d,%03d]_%d", det,ring,sec,strip,sample));
 135     TH1S* hSample = new TH1S(name.Data(),name.Data(), 1024,-.5,1023.5);
 136     hSample->SetXTitle("ADC");
 137     hSample->SetYTitle("Events");
 138     hSample->SetDirectory(0);
 139     sampleArray->AddAtAndExpand(hSample, sample);
 140   }
 141 }
 142
 143 //_____________________________________________________________________
 144 void AliFMDPedestalDA::AddSectorSummary(TObjArray* sectorArray,
 145                                         UShort_t det,
 146                                         Char_t   ring,
 147                                         UShort_t sec,
 148                                         UShort_t nStr)
 149 {
 150   TH1F* sumPed = new TH1F("Pedestals",
 151                           Form("Summary of pedestals in FMD%d%c[%02d]",
 152                                det, ring, sec),
 153                           nStr, -.5, nStr-.5);
 154   sumPed->SetXTitle("Strip");
 155   sumPed->SetYTitle("Pedestal [ADC]");
 156   sumPed->SetDirectory(0);
 157
 158   TH1F* sumNoise = static_cast<TH1F*>(sumPed->Clone("Noise"));
 159   sumNoise->SetYTitle("Noise [ADC]");
 160   sumNoise->SetDirectory(0);
 161
 162   Int_t n = sectorArray->GetEntriesFast();
 163   sectorArray->AddAtAndExpand(sumPed,   n + kPedestalOffset - 1);
 164   sectorArray->AddAtAndExpand(sumNoise, n + kNoiseOffset - 1);
 165 }
 166
 167 //_____________________________________________________________________
 168 void AliFMDPedestalDA::FillChannels(AliFMDDigit* digit)
 169 {
 170   // Fill ADC values from a digit into the corresponding histogram.
 171   //
 172   // Parameters:
 173   //    digit Digit to fill ADC values for.
 174   UShort_t det   = digit->Detector();
 175   Char_t   ring  = digit->Ring();
 176   UShort_t sec   = digit->Sector();
 177   UShort_t strip = digit->Strip();
 178
 179   AliFMDParameters* pars     = AliFMDParameters::Instance();
 180   UInt_t            samples  = pars->GetSampleRate(det, ring, sec, strip);
 181   for (UInt_t sample = 0; sample < samples; sample++) {
 182     TH1S* hSample = GetChannel(det, ring, sec, strip, sample);
 183     if (!hSample) continue;
 184
 185     hSample->Fill(digit->Count(sample));
 186   }
 187
 188 }
 189
 190 //_____________________________________________________________________
 191 void AliFMDPedestalDA::MakeSummary(UShort_t det, Char_t ring)
 192 {
 193   //Create summary hists for FMD pedestals
 194   // std::cout << "Making summary for FMD" << det << ring << " ..."
 195   //           << std::endl;
 196   switch (det) {
 197   case 1:
 198     fSummaryFMD1i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 199     break;
 200   case 2:
 201     switch (ring) {
 202     case 'I': case 'i':
 203       fSummaryFMD2i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 204       break;
 205     case 'O': case 'o':
 206       fSummaryFMD2o = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 207       break;
 208     }
 209     break;
 210   case 3:
 211     switch (ring) {
 212     case 'I': case 'i':
 213       fSummaryFMD3i = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 214       break;
 215     case 'O': case 'o':
 216       fSummaryFMD3o = MakeSummaryHistogram("ped", "Pedestals", det, ring);
 217       break;
 218     }
 219     break;
 220   }
 221 }
 222
 223 //_____________________________________________________________________
 224 void AliFMDPedestalDA::Analyse(UShort_t det,
 225                                Char_t   ring,
 226                                UShort_t sec,
 227                                UShort_t strip)
 228 {
 229   // Analyse a strip.  That is, compute the mean and spread of the ADC
 230   // spectra for all strips.  Also output on files the values.
 231   //
 232   // Parameters:
 233   //     det   Detector
 234   //     ring  Ring
 235   //     sec   Sector
 236   //     strip Strip.
 237   AliFMDParameters* pars     = AliFMDParameters::Instance();
 238   TH2* summary = 0;
 239   switch (det) {
 240   case 1: summary = fSummaryFMD1i; break;
 241   case 2:
 242     switch (ring) {
 243     case 'I':  summary = fSummaryFMD2i; break;
 244     case 'O':  summary = fSummaryFMD2o; break;
 245     }
 246     break;
 247   case 3:
 248     switch (ring) {
 249     case 'I':  summary = fSummaryFMD3i; break;
 250     case 'O':  summary = fSummaryFMD3o; break;
 251     }
 252     break;
 253   }
 254 #if 0
 255   static bool first = true;
 256   if (summary && first) {
 257     std::cout << "Filling summary " << summary->GetName() << std::endl;
 258     first = false;
 259   }
 260 #endif
 261   // Float_t           factor   = pars->GetPedestalFactor();
 262   UInt_t            samples  = pars->GetSampleRate(det, ring, sec, strip);
 263   for (UShort_t sample = 0; sample < samples; sample++) {
 264
 265     TH1S* hChannel = GetChannel(det, ring, sec, strip,sample);
 266     if(!hChannel || hChannel->GetEntries() == 0) {
 267       //AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
 268       //            det,ring,sec,strip));
 269       return;
 270     }
 271
 272     AliDebug(50, Form("Fitting FMD%d%c[%02d,%03d] with %d entries",
 273                       det,ring,sec,strip, int(hChannel->GetEntries())));
 274     TF1 fitFunc("fitFunc","gausn",0,300);
 275     fitFunc.SetParameters(100,100,1);
 276     hChannel->Fit("fitFunc","Q0+","",10,200);
 277
 278     Float_t mean = hChannel->GetMean();
 279     Float_t rms  = hChannel->GetRMS();
 280
 281
 282
 283     hChannel->GetXaxis()->SetRangeUser(mean-5*rms,mean+5*rms);
 284
 285     mean = hChannel->GetMean();
 286     rms  = hChannel->GetRMS();
 287
 288
 289     UShort_t ddl, board, altro, channel;
 290     UShort_t timebin;
 291
 292     pars->Detector2Hardware(det,ring,sec,strip,sample,
 293                             ddl,board,altro,channel,timebin);
 294     Int_t idx = HWIndex(ddl, board, altro, channel);
 295     if (idx >= 0) {
 296       fMinTimebin[idx] = TMath::Min(Short_t(timebin),   fMinTimebin[idx]);
 297       fMaxTimebin[idx] = TMath::Max(Short_t(timebin+1), fMaxTimebin[idx]);
 298     }
 299
 300     std::ostream* zsFile = 0;
 301     switch(det) {
 302     case 1:  zsFile = &fZSfileFMD1; break;
 303     case 2:  zsFile = &fZSfileFMD2; break;
 304     case 3:  zsFile = &fZSfileFMD3; break;
 305     default: AliWarning("Unknown sample!"); break;
 306
 307     }
 308     *zsFile << board   << ','
 309             << altro   << ','
 310             << channel << ','
 311             << timebin << ','
 312             << mean    << ','
 313             << rms     << "\n";
 314
 315     Float_t chi2ndf = 0;
 316
 317
 318     if(fitFunc.GetNDF())
 319       chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
 320
 321
 322     Int_t sampleToWrite = 2;
 323     if      (samples == 2) sampleToWrite = 1;
 324     else if (samples <  2) sampleToWrite = 0;
 325
 326     hChannel->GetXaxis()->SetRange(1,1024);
 327
 328     if(sample != sampleToWrite) continue;
 329
 330
 331     fOutputFile << det                         << ','
 332                 << ring                        << ','
 333                 << sec                         << ','
 334                 << strip                       << ','
 335                 << mean                        << ','
 336                 << rms                         << ','
 337                 << fitFunc.GetParameter(1)     << ','
 338                 << fitFunc.GetParameter(2)     << ','
 339                 << chi2ndf                     <<"\n";
 340
 341     if (summary) {
 342       Int_t bin = summary->FindBin(sec, strip);
 343       summary->SetBinContent(bin, mean);
 344       summary->SetBinError(bin, rms);
 345     }
 346
 347     if(fSaveHistograms  ) {
 348       TH1F* sumPed   = GetSectorSummary(det, ring, sec, true);
 349       TH1F* sumNoise = GetSectorSummary(det, ring, sec, false);
 350 #if 0
 351       gDirectory->cd(GetSectorPath(det, ring, sec, kTRUE));
 352       TH1F* sumPed   = dynamic_cast<TH1F*>(gDirectory->Get("Pedestals"));
 353       TH1F* sumNoise = dynamic_cast<TH1F*>(gDirectory->Get("Noise"));
 354       Int_t nStr = (ring == 'I' ? 512 : 256);
 355       if (!sumPed) {
 356         sumPed = new TH1F("Pedestals",
 357                           Form("Summary of pedestals in FMD%d%c[%02d]",
 358                                det, ring, sec),
 359                           nStr, -.5, nStr-.5);
 360         sumPed->SetXTitle("Strip");
 361         sumPed->SetYTitle("Pedestal [ADC]");
 362         sumPed->SetDirectory(gDirectory);
 363       }
 364       if (!sumNoise) {
 365         sumNoise = new TH1F("Noise",
 366                             Form("Summary of noise in FMD%d%c[%02d]",
 367                                  det, ring, sec),
 368                             nStr, -.5, nStr-.5);
 369         sumNoise->SetXTitle("Strip");
 370         sumNoise->SetYTitle("Noise [ADC]");
 371
 372         sumNoise->SetDirectory(gDirectory);
 373       }
 374 #endif
 375       sumPed->SetBinContent(strip+1, mean);
 376       sumPed->SetBinError(strip+1, rms);
 377       sumNoise->SetBinContent(strip+1, rms);
 378
 379 #if 0
 380       if(sumNoise->GetEntries() == nStr)
 381         sumNoise->Write(sumNoise->GetName(),TObject::kOverwrite);
 382       if(sumPed->GetEntries() == nStr)
 383         sumPed->Write(sumPed->GetName(),TObject::kOverwrite);
 384 #endif
 385
 386       fPedSummary.SetBinContent(fCurrentChannel,mean);
 387
 388       fNoiseSummary.SetBinContent(fCurrentChannel,rms);
 389       fCurrentChannel++;
 390
 391 #if 0
 392       gDirectory->cd(GetStripPath(det, ring, sec, strip, kTRUE));
 393       hChannel->GetXaxis()->SetRange(1,1024);
 394
 395       hChannel->Write();
 396 #endif
 397     }
 398   }
 399 }
 400
 401 //_____________________________________________________________________
 402 void AliFMDPedestalDA::Terminate(TFile* diagFile)
 403 {
 404   // Called at the end of a job.  Fills in missing time-bins and
 405   // closes output files
 406   if(fSaveHistograms && diagFile) {
 407     diagFile->cd();
 408
 409     fPedSummary.Write();
 410     fNoiseSummary.Write();
 411   }
 412   AliFMDAltroMapping* map = AliFMDParameters::Instance()->GetAltroMap();
 413   for (Int_t i = 0; i < 3; i++) {
 414     std::ofstream& out = (i == 0 ? fZSfileFMD1 :
 415                           i == 1 ? fZSfileFMD2 :
 416                           fZSfileFMD3);
 417     if (out.is_open() && fSeenDetectors[i]) {
 418       FillinTimebins(out, map->Detector2DDL(i+1));
 419     }
 420     if (!fSeenDetectors[i]) {
 421       TString n(Form("ddl%d.csv",3072+map->Detector2DDL(i+1)));
 422       gSystem->Unlink(n.Data());
 423     }
 424   }
 425
 426 }
 427
 428 //_____________________________________________________________________
 429 void AliFMDPedestalDA::FillinTimebins(std::ofstream& out, UShort_t /*ddl*/)
 430 {
 431   //
 432   // Fill missing timebins
 433   //
 434 #if 0
 435   unsigned short  boards[] = { 0x0, 0x1, 0x10, 0x11, 0xFFFF };
 436   unsigned short* board    = boards;
 437   while ((*boards) != 0xFFFF) {
 438     for (UShort_t altro = 0; altro < 3; altro++) {
 439       for (UShort_t channel = 0; channel < 16; channel++) {
 440         Int_t idx = HWIndex(ddl, *board, altro, channel);
 441         if (idx < 0) {
 442           AliWarning(Form("Invalid index for %4d/0x%02x/0x%x/0x%x: %d",
 443                           ddl, *board, altro, channel, idx));
 444           continue;
 445         }
 446         Short_t min = fMinTimebin[idx];
 447         Short_t max = fMaxTimebin[idx];
 448
 449         // Channel not seen at all.
 450         if (min > 1023 || max < 0) continue;
 451
 452         out << "# Extra timebins for 0x" << std::hex
 453             << board << ',' << altro << ',' << channel
 454             << " got time-bins " << min << " to " << max-1
 455             << std::dec << std::endl;
 456
 457         for (UShort_t t = 15; t < min; t++)
 458           // Write a phony line
 459           out << board << "," << altro << "," << channel << ","
 460               << t << "," << 1023 << "," << 0 << std::endl;
 461
 462         for (UShort_t t = max; t < 1024; t++)
 463           // Write a phony line
 464           out << board << "," << altro << "," << channel << ","
 465               << t << "," << 1023 << "," << 0 << std::endl;
 466       } // channel loop
 467     } // altro loop
 468   } // board loop
 469   // Write trailer, and close
 470 #endif
 471   out.write("# EOF\n", 6);
 472   out.close();
 473 }
 474
 475 //_____________________________________________________________________
 476 void AliFMDPedestalDA::WriteHeaderToFile()
 477 {
 478   //
 479   // Write headers to output files
 480   //
 481   AliFMDParameters* pars       = AliFMDParameters::Instance();
 482   fOutputFile.write(Form("# %s \n",pars->GetPedestalShuttleID()),13);
 483   TDatime now;
 484   fOutputFile << "# This file created from run # " << fRunno
 485               << " @ " << now.AsString() << std::endl;
 486   fOutputFile.write("# Detector, "
 487                     "Ring, "
 488                     "Sector, "
 489                     "Strip, "
 490                     "Pedestal, "
 491                     "Noise, "
 492                     "Mu, "
 493                     "Sigma, "
 494                     "Chi2/NDF \n", 71);
 495
 496   std::ostream* zss[] = { &fZSfileFMD1, &fZSfileFMD2, &fZSfileFMD3, 0 };
 497   for (size_t i = 0; i < 3; i++)  {
 498     *(zss[i]) << "# FMD " << (i+1) << " pedestals \n"
 499               << "# board, "
 500               << "altro, "
 501               << "channel, "
 502               << "timebin, "
 503               << "pedestal, "
 504               << "noise\n";
 505     *(zss[i]) << "# This file created from run # " << fRunno
 506               << " @ " << now.AsString() << std::endl;
 507   }
 508 }
 509
 510 //_____________________________________________________________________
 511 TH1S* AliFMDPedestalDA::GetChannel(UShort_t det,
 512                                    Char_t   ring,
 513                                    UShort_t sec,
 514                                    UShort_t strip,
 515                                    UInt_t   sample)
 516 {
 517   // Get the histogram corresponding to a strip sample.
 518   //
 519   // Parameters:
 520   //     det    Detector
 521   //     ring   Ring
 522   //     sec    Sector
 523   //     strip  Strip
 524   //     sample Sample
 525   //
 526   // Return:
 527   //     ADC spectra of a strip.
 528   TObjArray* sampleArray = GetStripArray(det, ring, sec, strip);
 529   if (!sampleArray) return 0;
 530   TH1S*      hSample     = static_cast<TH1S*>(sampleArray->At(sample));
 531   if (!hSample) {
 532     AliErrorF("No channel histogram for FMD%d%c[%02d,%03d]_%d",
 533               det, ring, sec, strip, sample);
 534     sampleArray->ls();
 535     AliErrorF("Path is %s <- %s <- %s <- %s",
 536               sampleArray->GetName(),
 537               GetSectorArray(det, ring, sec)->GetName(),
 538               GetRingArray(det, ring)->GetName(),
 539               GetDetectorArray(det)->GetName());
 540
 541   }
 542   return hSample;
 543
 544 }
 545 //_____________________________________________________________________
 546 TH1F* AliFMDPedestalDA::GetSectorSummary(UShort_t det,
 547                                          Char_t   ring,
 548                                          UShort_t sec,
 549                                          Bool_t   pedNotNoise)
 550 {
 551   TObjArray* secArray    = GetSectorArray(det, ring, sec);
 552   Int_t      n           = secArray->GetEntriesFast();
 553   Int_t      i           = n - (pedNotNoise ? kNoiseOffset : kPedestalOffset);
 554   return static_cast<TH1F*>(secArray->At(i));
 555 }
 556
 557 //_____________________________________________________________________
 558 //
 559 //EOF
 560 //