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