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 "TH1.h"
  37 #include "TMath.h"
  38 #include "TGraphErrors.h"
  39 #include "AliFMDParameters.h"
  40
  41 //_____________________________________________________________________
  42 ClassImp(AliFMDGainDA)
  43 #if 0 // Do not delete - here to let Emacs indent properly
  44 ;
  45 #endif
  46
  47 //_____________________________________________________________________
  48 AliFMDGainDA::AliFMDGainDA()
  49   : AliFMDBaseDA(),
  50     fGainArray(),
  51     fHighPulse(256),
  52     fEventsPerChannel(10),
  53     fCurrentPulse(10),
  54     fCurrentChannel(10),
  55     fNumberOfStripsPerChip(128),
  56     fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
  57     fCurrentSummaryStrip(1)
  58 {
  59   fCurrentPulse.Reset(0);
  60   fCurrentChannel.Reset(0);
  61   fOutputFile.open("gains.csv");
  62   fGainArray.SetOwner();
  63 }
  64
  65 //_____________________________________________________________________
  66 AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA)
  67   :  AliFMDBaseDA(gainDA),
  68      fGainArray(gainDA.fGainArray),
  69      fHighPulse(gainDA.fHighPulse),
  70      fEventsPerChannel(gainDA.fEventsPerChannel),
  71      fCurrentPulse(gainDA.fCurrentPulse),
  72      fCurrentChannel(gainDA.fCurrentChannel),
  73      fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
  74      fSummaryGains(gainDA.fSummaryGains),
  75      fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip)
  76 {
  77   fCurrentPulse.Reset(0);
  78   fCurrentChannel.Reset(0);
  79 }
  80
  81 //_____________________________________________________________________
  82 AliFMDGainDA::~AliFMDGainDA()
  83 {
  84 }
  85
  86 //_____________________________________________________________________
  87 void AliFMDGainDA::Init()
  88 {
  89
  90
  91   Int_t nEventsRequired = 0;
  92
  93   //for(UShort_t det=1; det<=3;det++) {
  94   //  UShort_t firstring = (det == 1 ? 1 : 0);
  95   //  for(UShort_t iring = firstring; iring <=1;iring++) {
  96   //    Char_t ring = (iring == 1 ? 'I' : 'O');
  97   //    for(UShort_t board =0 ; board <=1; board++) {
  98   //    Int_t idx = GetHalfringIndex(det,ring,board);
  99   for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++)
 100     {
 101
 102       Int_t nEvents = 0;
 103       if(fPulseSize.At(idx))
 104         nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
 105       fEventsPerChannel.AddAt(nEvents,idx);
 106       if(nEvents>nEventsRequired) nEventsRequired = nEvents * fNumberOfStripsPerChip;
 107
 108     }
 109   //}
 110   // }
 111
 112   //8 pulser values * 128 strips * 100 samples
 113
 114
 115   SetRequiredEvents(nEventsRequired);
 116
 117   TObjArray* detArray;
 118   TObjArray* ringArray;
 119   TObjArray* sectorArray;
 120
 121   for(UShort_t det=1;det<=3;det++) {
 122     detArray = new TObjArray();
 123     detArray->SetOwner();
 124     fGainArray.AddAtAndExpand(detArray,det);
 125     for (UShort_t ir = 0; ir < 2; ir++) {
 126       Char_t   ring = (ir == 0 ? 'O' : 'I');
 127       UShort_t nsec = (ir == 0 ? 40  : 20);
 128       UShort_t nstr = (ir == 0 ? 2 : 4);
 129       ringArray = new TObjArray();
 130       ringArray->SetOwner();
 131       detArray->AddAtAndExpand(ringArray,ir);
 132       for(UShort_t sec =0; sec < nsec;  sec++)  {
 133         sectorArray = new TObjArray();
 134         sectorArray->SetOwner();
 135         ringArray->AddAtAndExpand(sectorArray,sec);
 136         for(UShort_t strip = 0; strip < nstr; strip++) {
 137           TH1S* hChannel = new TH1S(Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
 138                                     Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
 139                                     1024,0,1023);
 140           hChannel->SetDirectory(0);
 141           sectorArray->AddAtAndExpand(hChannel,strip);
 142         }
 143       }
 144     }
 145   }
 146 }
 147
 148 //_____________________________________________________________________
 149 void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray,
 150                                        UShort_t det  ,
 151                                        Char_t   ring,
 152                                        UShort_t sec,
 153                                        UShort_t strip)
 154 {
 155   TGraphErrors* hChannel  = new TGraphErrors();
 156   hChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
 157   hChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC",
 158                           det, ring, sec, strip));
 159   sectorArray->AddAtAndExpand(hChannel,strip);
 160 }
 161
 162 //_____________________________________________________________________
 163 void AliFMDGainDA::FillChannels(AliFMDDigit* digit) {
 164
 165   UShort_t det   = digit->Detector();
 166   Char_t   ring  = digit->Ring();
 167   UShort_t sec   = digit->Sector();
 168   UShort_t strip = digit->Strip();
 169
 170   //Strip is always seen as the first in a VA chip. All other strips are junk
 171   if(strip % fNumberOfStripsPerChip) return;
 172   Int_t vaChip   = strip / fNumberOfStripsPerChip;
 173   TH1S* hChannel = GetChannelHistogram(det, ring, sec, vaChip);
 174
 175   hChannel->Fill(digit->Counts());
 176   UpdatePulseAndADC(det,ring,sec,strip);
 177 }
 178
 179 //_____________________________________________________________________
 180 void AliFMDGainDA::Analyse(UShort_t det,
 181                            Char_t   ring,
 182                            UShort_t sec,
 183                            UShort_t strip) {
 184   TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
 185   if(!grChannel->GetN()) {
 186     // AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
 187     //                 det, ring , sec, strip));
 188     return;
 189   }
 190   TF1 fitFunc("fitFunc","pol1",-10,280);
 191   fitFunc.SetParameters(100,3);
 192   grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
 193
 194   Float_t gain    = -1;
 195   Float_t error   = -1;
 196   Float_t chi2ndf = -1;
 197   if((fitFunc.GetParameter(1)) == (fitFunc.GetParameter(1))) {
 198     gain    = fitFunc.GetParameter(1);
 199     error   = fitFunc.GetParError(1);
 200     if(fitFunc.GetNDF())
 201       chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
 202   }
 203
 204   fOutputFile << det                         << ','
 205               << ring                        << ','
 206               << sec                         << ','
 207               << strip                       << ','
 208               << gain                        << ','
 209               << error                       << ','
 210               << chi2ndf                     <<"\n";
 211
 212   //due to RCU trouble, first strips on VAs are excluded
 213   if(strip%128 != 0) {
 214
 215     fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
 216     fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
 217
 218     fCurrentSummaryStrip++;
 219   }
 220   if(fSaveHistograms) {
 221     gDirectory->cd(GetSectorPath(det,ring, sec, kTRUE));
 222
 223     TH1F* summary = dynamic_cast<TH1F*>(gDirectory->Get("Summary"));
 224     if (!summary) {
 225       Int_t nStr = (ring == 'I' ? 512 : 256);
 226       summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]",
 227                                          det, ring, sec),
 228                          nStr, -.5, nStr-.5);
 229       summary->SetXTitle("Strip");
 230       summary->SetYTitle("Gain [ADC/DAC]");
 231       summary->SetDirectory(gDirectory);
 232     }
 233     summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
 234     summary->SetBinError(strip+1, fitFunc.GetParError(1));
 235
 236     gDirectory->cd(GetStripPath(det,ring,sec,strip, kTRUE));
 237     grChannel->SetName(Form("FMD%d%c[%02d,%03d]",det,ring,sec,strip));
 238     // grChannel->SetDirectory(gDirectory);
 239     grChannel->Write();
 240     // grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
 241   }
 242 }
 243
 244 //_____________________________________________________________________
 245 void AliFMDGainDA::Terminate(TFile* diagFile)
 246 {
 247   diagFile->cd();
 248   fSummaryGains.Write();
 249 }
 250
 251 //_____________________________________________________________________
 252 void AliFMDGainDA::WriteHeaderToFile()
 253 {
 254   AliFMDParameters* pars       = AliFMDParameters::Instance();
 255   fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
 256   fOutputFile.write("# Detector, "
 257                     "Ring, "
 258                     "Sector, "
 259                     "Strip, "
 260                     "Gain, "
 261                     "Error, "
 262                     "Chi2/NDF \n",56);
 263
 264 }
 265
 266 //_____________________________________________________________________
 267 TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
 268                                         Char_t   ring,
 269                                         UShort_t sec,
 270                                         UShort_t strip)
 271 {
 272
 273   UShort_t  Ring = 1;
 274   if(ring == 'O')
 275     Ring = 0;
 276
 277
 278   TObjArray* detArray  = static_cast<TObjArray*>(fGainArray.At(det));
 279   TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
 280   TObjArray* secArray  = static_cast<TObjArray*>(ringArray->At(sec));
 281   TH1S* hChannel       = static_cast<TH1S*>(secArray->At(strip));
 282
 283   return hChannel;
 284 }
 285
 286 //_____________________________________________________________________
 287 TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
 288                                        Char_t   ring,
 289                                        UShort_t sec,
 290                                        UShort_t strip)
 291 {
 292   UShort_t      iring     = (ring == 'O' ? 0 : 1);
 293   TObjArray*    detArray  = static_cast<TObjArray*>(fDetectorArray.At(det));
 294   TObjArray*    ringArray = static_cast<TObjArray*>(detArray->At(iring));
 295   TObjArray*    secArray  = static_cast<TObjArray*>(ringArray->At(sec));
 296   TGraphErrors* hChannel  = static_cast<TGraphErrors*>(secArray->At(strip));
 297
 298   return hChannel;
 299 }
 300
 301 //_____________________________________________________________________
 302 void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
 303                                      Char_t ring,
 304                                      UShort_t sec,
 305                                      UShort_t strip)
 306 {
 307
 308   AliFMDParameters* pars = AliFMDParameters::Instance();
 309   UInt_t ddl, board,chip,ch;
 310   pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,ch);
 311   Int_t halfring = GetHalfringIndex(det,ring,board/16);
 312
 313   if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
 314     return;
 315   if(strip % fNumberOfStripsPerChip) return;
 316   if(((GetCurrentEvent()) % fPulseLength.At(halfring)) && GetCurrentEvent() > 0) return;
 317
 318   Int_t vaChip = strip/fNumberOfStripsPerChip;
 319   TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
 320
 321   if(!hChannel->GetEntries()) {
 322     AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
 323                     det, ring , sec, strip));
 324     return;
 325   }
 326   Double_t mean      = hChannel->GetMean();
 327   Double_t rms       = hChannel->GetRMS();
 328   Double_t pulse     = Double_t(fCurrentPulse.At(halfring)) * fPulseSize.At(halfring);
 329   Int_t    firstBin  = hChannel->GetXaxis()->GetFirst();
 330   Int_t    lastBin   = hChannel->GetXaxis()->GetLast();
 331   hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
 332
 333   mean               = hChannel->GetMean();
 334   rms                = hChannel->GetRMS();
 335
 336   hChannel->GetXaxis()->SetRange(firstBin,lastBin);
 337
 338   Int_t channelNumber = strip + (GetCurrentEvent()-1)/((fPulseLength.At(halfring)*fHighPulse)/fPulseSize.At(halfring));
 339
 340   TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
 341
 342   channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
 343   channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
 344
 345   if(fSaveHistograms) {
 346     gDirectory->cd(GetStripPath(det,ring,sec,channelNumber));
 347     hChannel->Write(Form("%s_pulse_%03d",hChannel->GetName(),(Int_t)pulse));
 348
 349   }
 350
 351   hChannel->Reset();
 352
 353 }
 354
 355 //_____________________________________________________________________
 356 void AliFMDGainDA::ResetPulseAndUpdateChannel()
 357 {
 358   fCurrentPulse.Reset(0);
 359 }
 360
 361 //_____________________________________________________________________
 362 void AliFMDGainDA::FinishEvent()
 363 {
 364   for(UShort_t det=1; det<=3;det++) {
 365     UShort_t firstring = (det == 1 ? 1 : 0);
 366     for(UShort_t iring = firstring; iring <=1;iring++) {
 367       Char_t ring = (iring == 1 ? 'I' : 'O');
 368       for(UShort_t board =0 ; board <=1; board++) {
 369         Int_t idx = GetHalfringIndex(det,ring,board);
 370
 371         if( !fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
 372           continue;
 373         if(GetCurrentEvent()>0 && ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
 374           fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
 375
 376         if(GetCurrentEvent()>0 && ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
 377           fCurrentPulse.AddAt(0,idx);
 378       }
 379     }
 380   }
 381 }
 382 //_____________________________________________________________________
 383 //
 384 //EOF
 385 //