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