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