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