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