FMD/AliFMDGainDA.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    AliFMDGainDA.cxx
  17     @author  Hans Hjersing Dalsgaard <canute@nbi.dk>
  18     @date    Mon Mar 13 13:46:05 2008
  19     @brief   Derived class for the pulse gain detector algorithm.
  20 */
  21 //____________________________________________________________________
  22 //
  23 // This class contains the implementation of the gain detector
  24 // algorithms (DA) for the FMD.  The data is collected in histograms
  25 // that are reset for each pulse length after the mean and standard
  26 // deviation are put into a TGraphErrors object. After a certain
  27 // number of pulses (usually 8) the graph is fitted to a straight
  28 // line. The gain is then slope of this line as it combines the known
  29 // pulse and the response of the detector.
  30 //
  31 #include "AliFMDGainDA.h"
  32 #include "AliFMDAltroMapping.h"
  33 #include "AliFMDParameters.h"
  34 #include "AliFMDCalibGain.h"
  35 #include "AliFMDDigit.h"
  36 #include "AliLog.h"
  37 #include <TFile.h>
  38 #include <TF1.h>
  39 #include <TH1S.h>
  40 #include <TGraphErrors.h>
  41 #include <TDatime.h>
  42 #include <TH2.h>
  43
  44 //_____________________________________________________________________
  45 ClassImp(AliFMDGainDA)
  46 #if 0 // Do not delete - here to let Emacs indent properly
  47 ;
  48 #endif
  49
  50 //_____________________________________________________________________
  51 AliFMDGainDA::AliFMDGainDA()
  52   : AliFMDBaseDA(),
  53     fHighPulse(256),
  54     fEventsPerChannel(10),
  55     fCurrentPulse(10),
  56     fCurrentChannel(10),
  57     fNumberOfStripsPerChip(128),
  58     fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
  59     fCurrentSummaryStrip(1),
  60     fGainFMD1i(0),
  61     fGainFMD2i(0),
  62     fGainFMD2o(0),
  63     fGainFMD3i(0),
  64     fGainFMD3o(0),
  65     fChi2FMD1i(0),
  66     fChi2FMD2i(0),
  67     fChi2FMD2o(0),
  68     fChi2FMD3i(0),
  69     fChi2FMD3o(0)
  70 {
  71   // Constructor
  72   //
  73   // Parameters:
  74   //   None
  75   fCurrentPulse.Reset(0);
  76   fCurrentChannel.Reset(0);
  77   fOutputFile.open("gains.csv");
  78   fDiagnosticsFilename = "diagnosticsGain.root";
  79 }
  80
  81 //_____________________________________________________________________
  82 AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA)
  83   : AliFMDBaseDA(gainDA),
  84     fHighPulse(gainDA.fHighPulse),
  85     fEventsPerChannel(gainDA.fEventsPerChannel),
  86     fCurrentPulse(gainDA.fCurrentPulse),
  87     fCurrentChannel(gainDA.fCurrentChannel),
  88     fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
  89     fSummaryGains(gainDA.fSummaryGains),
  90     fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip),
  91     fGainFMD1i(gainDA.fGainFMD1i),
  92     fGainFMD2i(gainDA.fGainFMD2i),
  93     fGainFMD2o(gainDA.fGainFMD2o),
  94     fGainFMD3i(gainDA.fGainFMD3i),
  95     fGainFMD3o(gainDA.fGainFMD3o),
  96     fChi2FMD1i(gainDA.fChi2FMD1i),
  97     fChi2FMD2i(gainDA.fChi2FMD2i),
  98     fChi2FMD2o(gainDA.fChi2FMD2o),
  99     fChi2FMD3i(gainDA.fChi2FMD3i),
 100     fChi2FMD3o(gainDA.fChi2FMD3o)
 101 {
 102   // Copy Constructor
 103   //
 104   // Parameters:
 105   //   gainDA   Object to copy from
 106   fCurrentPulse.Reset(0);
 107   fCurrentChannel.Reset(0);
 108 }
 109
 110 //_____________________________________________________________________
 111 AliFMDGainDA::~AliFMDGainDA()
 112 {
 113   // Destructor
 114   //
 115   // Parameters:
 116   //   None
 117 }
 118
 119 //_____________________________________________________________________
 120 void AliFMDGainDA::Init()
 121 {
 122   // Initialize
 123   //
 124   // Parameters:
 125   //   None
 126   Int_t nEventsRequired = 0;
 127
 128   for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++) {
 129     Int_t nEvents = 0;
 130     if(fPulseSize.At(idx))
 131       nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
 132     fEventsPerChannel.AddAt(nEvents,idx);
 133     if(nEvents>nEventsRequired)
 134       nEventsRequired = nEvents * fNumberOfStripsPerChip;
 135   }
 136   SetRequiredEvents(nEventsRequired);
 137 }
 138
 139 //_____________________________________________________________________
 140 void AliFMDGainDA::InitContainer(TDirectory* dir)
 141 {
 142   AliFMDBaseDA::InitContainer(dir);
 143
 144   for(UShort_t det=1;det<=3;det++) {
 145     UShort_t sr = (det == 1 ? 1 : 0);
 146     for (UShort_t ir = sr; ir < 2; ir++) {
 147       Char_t   ring = (ir == 0 ? 'O' : 'I');
 148       UShort_t nsec = (ir == 0 ? 40  : 20);
 149       UShort_t nva  = (ir == 0 ? 2 : 4);
 150       for(UShort_t sec =0; sec < nsec;  sec++)  {
 151         TObjArray* sectorArray = GetSectorArray(det, ring, sec);
 152         TObjArray* cache = new TObjArray(nva);
 153         cache->SetName("Cache");
 154         cache->SetOwner();
 155         Int_t n = sectorArray->GetEntriesFast();
 156         sectorArray->AddAtAndExpand(cache, n);
 157         for(UShort_t va = 0; va < nva; va++) {
 158           TH1S* hChannel = new TH1S(Form("FMD%d%c[%02d]_va%d",
 159                                          det,ring,sec,va),
 160                                     Form("FMD%d%c[%02d] VA%d Cache",
 161                                          det,ring,sec,va),
 162                                     1024,-.5,1023.5);
 163           hChannel->SetDirectory(0);
 164           cache->AddAtAndExpand(hChannel,va);
 165         }
 166       }
 167     }
 168   }
 169 }
 170
 171 //_____________________________________________________________________
 172 void AliFMDGainDA::AddChannelContainer(TObjArray* stripArray,
 173                                        UShort_t det  ,
 174                                        Char_t   ring,
 175                                        UShort_t sec,
 176                                        UShort_t strip)
 177 {
 178   // Make a channel container
 179   //
 180   // Parameters:
 181   //  sectorArray    Sectors
 182   //  det            Detector number
 183   //  ring           Ring identifier
 184   //  sec            Sector number
 185   //  strip          Strip number
 186
 187   AliFMDParameters* pars     = AliFMDParameters::Instance();
 188   UShort_t          board    = pars->GetAltroMap()->Sector2Board(ring, sec);
 189   Int_t             halfring = GetHalfringIndex(det,ring,board/16);
 190   Int_t             dPulse   = fPulseSize.At(halfring);
 191   Int_t             nPulses  = dPulse > 0 ? fHighPulse / dPulse : 0;
 192
 193   TGraphErrors* gChannel  = new TGraphErrors(nPulses);
 194   gChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
 195   gChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC",
 196                           det, ring, sec, strip));
 197   stripArray->AddAtAndExpand(gChannel,0);
 198 }
 199
 200 //_____________________________________________________________________
 201 void AliFMDGainDA::AddSectorSummary(TObjArray* sectorArray,
 202                                     UShort_t det,
 203                                     Char_t   ring,
 204                                     UShort_t sec,
 205                                     UShort_t nStr)
 206 {
 207   TH1F* summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]",
 208                                            det, ring, sec),
 209                            nStr, -.5, nStr-.5);
 210   summary->SetXTitle("Strip");
 211   summary->SetYTitle("Gain [ADC/DAC]");
 212   summary->SetDirectory(0);
 213
 214   Int_t n = sectorArray->GetEntriesFast();
 215   sectorArray->AddAtAndExpand(summary, n);
 216 }
 217
 218 //_____________________________________________________________________
 219 void AliFMDGainDA::FillChannels(AliFMDDigit* digit)
 220 {
 221   // Fill data into histogram
 222   //
 223   // Parameters:
 224   //   digit  Digit to get the data from
 225
 226   UShort_t det   = digit->Detector();
 227   Char_t   ring  = digit->Ring();
 228   UShort_t sec   = digit->Sector();
 229   UShort_t strip = digit->Strip();
 230
 231   //Strip is always seen as the first in a VA chip. All other strips are junk.
 232   //Strips are counted from zero on even sectors and from 511 on odd sectors...
 233
 234   if((sec%2)     && ((strip+1) % fNumberOfStripsPerChip)) return;
 235   if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
 236
 237   UShort_t vaChip   = strip / fNumberOfStripsPerChip;
 238   TH1S*    hChannel = GetChannelHistogram(det, ring, sec, vaChip);
 239   hChannel->Fill(digit->Counts());
 240   UpdatePulseAndADC(det,ring,sec,strip);
 241 }
 242
 243 //_____________________________________________________________________
 244 void AliFMDGainDA::MakeSummary(UShort_t det, Char_t ring)
 245 {
 246   //
 247   // Create summary hists for FMD gains and chi2 of the fits
 248   //
 249   switch (det) {
 250   case 1:
 251     fGainFMD1i = MakeSummaryHistogram("gain", "Gains", det, ring);
 252     fChi2FMD1i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
 253     break;
 254   case 2:
 255     switch (ring) {
 256     case 'I': case 'i':
 257       fGainFMD2i = MakeSummaryHistogram("gain", "Gains", det, ring);
 258       fChi2FMD2i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
 259       break;
 260     case 'O': case 'o':
 261       fGainFMD2o = MakeSummaryHistogram("gain", "Gains", det, ring);
 262       fChi2FMD2o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
 263       break;
 264     }
 265     break;
 266   case 3:
 267     switch (ring) {
 268     case 'I': case 'i':
 269       fGainFMD3i = MakeSummaryHistogram("gain", "Gains", det, ring);
 270       fChi2FMD3i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
 271       break;
 272     case 'O': case 'o':
 273       fGainFMD3o = MakeSummaryHistogram("gain", "Gains", det, ring);
 274       fChi2FMD3o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
 275       break;
 276     }
 277     break;
 278   }
 279 }
 280
 281 //_____________________________________________________________________
 282 void AliFMDGainDA::Analyse(UShort_t det,
 283                            Char_t   ring,
 284                            UShort_t sec,
 285                            UShort_t strip)
 286 {
 287   // Analyse result of a single strip
 288   //
 289   // Parameters:
 290   //  det            Detector number
 291   //  ring           Ring identifier
 292   //  sec            Sector number
 293   //  strip          Strip number
 294   TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
 295   if(!grChannel->GetN()) {
 296     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
 297                      det, ring , sec, strip));
 298     return;
 299   }
 300   TF1 fitFunc("fitFunc","pol1",-10,280);
 301   fitFunc.SetParameters(100,3);
 302   grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
 303
 304   Float_t gain    = -1;
 305   Float_t error   = -1;
 306   Float_t chi2ndf = -1;
 307   if((fitFunc.GetParameter(1)) == (fitFunc.GetParameter(1))) {
 308     gain    = fitFunc.GetParameter(1);
 309     error   = fitFunc.GetParError(1);
 310     if(fitFunc.GetNDF())
 311       chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
 312   }
 313
 314   fOutputFile << det                         << ','
 315               << ring                        << ','
 316               << sec                         << ','
 317               << strip                       << ','
 318               << gain                        << ','
 319               << error                       << ','
 320               << chi2ndf                     <<"\n";
 321
 322   //due to RCU trouble, first strips on VAs are excluded
 323   // if(strip%128 != 0) {
 324
 325   fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
 326   fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
 327
 328   fCurrentSummaryStrip++;
 329
 330   TH2* hGain = 0;
 331   TH2* hChi2 = 0;
 332   switch (det) {
 333   case 1: hGain = fGainFMD1i; hChi2 = fChi2FMD1i; break;
 334   case 2:
 335     switch (ring) {
 336     case 'I':  hGain = fGainFMD2i; hChi2 = fChi2FMD2i; break;
 337     case 'O':  hGain = fGainFMD2o; hChi2 = fChi2FMD2o; break;
 338     }
 339     break;
 340   case 3:
 341     switch (ring) {
 342     case 'I':  hGain = fGainFMD3i; hChi2 = fChi2FMD3i; break;
 343     case 'O':  hGain = fGainFMD3o; hChi2 = fChi2FMD3o; break;
 344     }
 345     break;
 346   }
 347   if (hGain && hChi2) {
 348     Int_t bin = hGain->FindBin(sec, strip);
 349     hGain->SetBinContent(bin, gain);
 350     hGain->SetBinError(bin, error);
 351     hChi2->SetBinContent(bin, chi2ndf);
 352   }
 353
 354     // }
 355   if(fSaveHistograms) {
 356     TH1F* summary = GetSectorSummary(det, ring, sec);
 357     summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
 358     summary->SetBinError(strip+1, fitFunc.GetParError(1));
 359   }
 360 }
 361
 362 //_____________________________________________________________________
 363 void AliFMDGainDA::Terminate(TFile* diagFile)
 364 {
 365   // End of file
 366   //
 367   // Parameters:
 368   //   None
 369   if(diagFile) {
 370     diagFile->cd();
 371     fSummaryGains.Write();
 372   }
 373 }
 374
 375 //_____________________________________________________________________
 376 void AliFMDGainDA::WriteHeaderToFile()
 377 {
 378   // Write header to the output file
 379   //
 380   // Parameters:
 381   //   None
 382   AliFMDParameters* pars       = AliFMDParameters::Instance();
 383   fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
 384   TDatime now;
 385   fOutputFile << "# This file created from run # " << fRunno
 386               << " @ " << now.AsString() << std::endl;
 387   fOutputFile.write("# Detector, "
 388                     "Ring, "
 389                     "Sector, "
 390                     "Strip, "
 391                     "Gain, "
 392                     "Error, "
 393                     "Chi2/NDF \n",56);
 394 }
 395
 396 //_____________________________________________________________________
 397 TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
 398                                         Char_t   ring,
 399                                         UShort_t sec,
 400                                         UShort_t va)
 401 {
 402   // Get the current histogram of a single strip
 403   //
 404   // Parameters:
 405   //  det            Detector number
 406   //  ring           Ring identifier
 407   //  sec            Sector number
 408   //  va             VA number
 409   TObjArray* secArray  = GetSectorArray(det, ring, sec);
 410   Int_t      n         = secArray->GetEntriesFast();
 411   TObjArray* cache     = static_cast<TObjArray*>(secArray->At(n-1));
 412   TH1S* hChannel       = static_cast<TH1S*>(cache->At(va));
 413
 414   return hChannel;
 415 }
 416
 417 //_____________________________________________________________________
 418 TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
 419                                        Char_t   ring,
 420                                        UShort_t sec,
 421                                        UShort_t strip)
 422 {
 423   // Get the graph of a single strip
 424   //
 425   // Parameters:
 426   //  det            Detector number
 427   //  ring           Ring identifier
 428   //  sec            Sector number
 429   //  strip          Strip number
 430   TObjArray*    stripArray = GetStripArray(det, ring, sec, strip);
 431   TGraphErrors* hChannel   = static_cast<TGraphErrors*>(stripArray->At(0));
 432
 433   return hChannel;
 434 }
 435
 436 //_____________________________________________________________________
 437 TH1F* AliFMDGainDA::GetSectorSummary(UShort_t det,
 438                                      Char_t   ring,
 439                                      UShort_t sec)
 440 {
 441   TObjArray* secArray    = GetSectorArray(det, ring, sec);
 442   Int_t      n           = secArray->GetEntriesFast();
 443   return static_cast<TH1F*>(secArray->At(n-2)); // Cache added later
 444 }
 445
 446 //_____________________________________________________________________
 447 void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
 448                                      Char_t ring,
 449                                      UShort_t sec,
 450                                      UShort_t strip)
 451 {
 452   // Go to next pulse size
 453   //
 454   // Parameters:
 455   //  det            Detector number
 456   //  ring           Ring identifier
 457   //  sec            Sector number
 458   //  strip          Strip number
 459
 460   AliFMDParameters* pars     = AliFMDParameters::Instance();
 461   UShort_t          board    = pars->GetAltroMap()->Sector2Board(ring, sec);
 462   Int_t             halfring = GetHalfringIndex(det,ring,board/16);
 463
 464   if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
 465     return;
 466
 467   if ((sec       % 2) && ((strip+1) % fNumberOfStripsPerChip)) return;
 468   if (((sec + 1) % 2) && (strip     % fNumberOfStripsPerChip)) return;
 469
 470   if(((GetCurrentEvent()) % fPulseLength.At(halfring))
 471      && GetCurrentEvent() > 0) return;
 472
 473   Int_t vaChip   = strip/fNumberOfStripsPerChip;
 474   TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
 475
 476   if(!hChannel->GetEntries()) {
 477     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
 478                     det, ring , sec, strip));
 479     return;
 480   }
 481   Double_t mean      = hChannel->GetMean();
 482   Double_t rms       = hChannel->GetRMS();
 483   Double_t pulse     = (Double_t(fCurrentPulse.At(halfring))
 484                         * fPulseSize.At(halfring));
 485   Int_t    firstBin  = hChannel->GetXaxis()->GetFirst();
 486   Int_t    lastBin   = hChannel->GetXaxis()->GetLast();
 487   hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
 488
 489   mean               = hChannel->GetMean();
 490   rms                = hChannel->GetRMS();
 491
 492   hChannel->GetXaxis()->SetRange(firstBin,lastBin);
 493
 494   Int_t         channelOff = ((GetCurrentEvent()-1)
 495                               / ((fPulseLength.At(halfring)*fHighPulse)
 496                                  / fPulseSize.At(halfring)));
 497   Int_t         channelNo  = (strip + ((sec % 2 == 1) ? -1 : 1) * channelOff);
 498   TGraphErrors* channel    = GetChannel(det,ring,sec,channelNo);
 499
 500   channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
 501   channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
 502
 503   if(fSaveHistograms) {
 504     TH1S* out =
 505       static_cast<TH1S*>(hChannel->Clone(Form("FMD%d%c[%02d,%03d]_0x%02x",
 506                                               det, ring, sec, channelNo,
 507                                               int(pulse))));
 508     out->SetTitle(Form("FMD%d%c[%02d,%03d] DAC=0x%02x (%3d)",
 509                        det, ring, sec, channelNo, int(pulse), int(pulse)));
 510     out->SetDirectory(0);
 511     TObjArray* arr = GetStripArray(det, ring, sec, channelNo);
 512     arr->AddAtAndExpand(out, fCurrentPulse.At(halfring)+1);
 513   }
 514
 515   hChannel->Reset();
 516
 517 }
 518
 519 //_____________________________________________________________________
 520 void AliFMDGainDA::ResetPulseAndUpdateChannel()
 521 {
 522   // Reset all
 523   //
 524   // Parameters:
 525   //   None
 526   fCurrentPulse.Reset(0);
 527 }
 528
 529 //_____________________________________________________________________
 530 void AliFMDGainDA::FinishEvent()
 531 {
 532   // End of event
 533   //
 534   // Parameters:
 535   //   None
 536   for(UShort_t det=1; det<=3;det++) {
 537     UShort_t firstring = (det == 1 ? 1 : 0);
 538     for(UShort_t iring = firstring; iring <=1;iring++) {
 539       Char_t ring = (iring == 1 ? 'I' : 'O');
 540       for(UShort_t board =0 ; board <=1; board++) {
 541         Int_t idx = GetHalfringIndex(det,ring,board);
 542
 543         if(!fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
 544           continue;
 545         if(GetCurrentEvent()>0 &&
 546            ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
 547           fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
 548
 549         if(GetCurrentEvent()>0 &&
 550            ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
 551           fCurrentPulse.AddAt(0,idx);
 552       }
 553     }
 554   }
 555 }
 556 //_____________________________________________________________________
 557 //
 558 //EOF
 559 //