Fixes for FMD DAs
[u/mrichter/AliRoot.git] / FMD / FMDutil / AliFMDGainDA.cxx
CommitLineData
f631f26d 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*/
cf291b42 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//
f631f26d 31#include "AliFMDGainDA.h"
b995fc28 32#include "AliFMDAltroMapping.h"
09b6c804 33#include "AliFMDParameters.h"
34#include "AliFMDCalibGain.h"
35#include "AliFMDDigit.h"
36#include "AliLog.h"
37#include <TFile.h>
38#include <TF1.h>
39#include <TH1S.h>
40#include <TGraphErrors.h>
6cf6e7a0 41#include <TDatime.h>
2a082c96 42#include <TH2.h>
f8022830 43#include <TROOT.h>
f631f26d 44
45//_____________________________________________________________________
46ClassImp(AliFMDGainDA)
cf291b42 47#if 0 // Do not delete - here to let Emacs indent properly
48;
49#endif
f631f26d 50
51//_____________________________________________________________________
cf291b42 52AliFMDGainDA::AliFMDGainDA()
53 : AliFMDBaseDA(),
cf291b42 54 fHighPulse(256),
22017a7e 55 fEventsPerChannel(10),
56 fCurrentPulse(10),
57 fCurrentChannel(10),
3bae5d02 58 fNumberOfStripsPerChip(128),
59 fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
2a082c96 60 fCurrentSummaryStrip(1),
61 fGainFMD1i(0),
62 fGainFMD2i(0),
63 fGainFMD2o(0),
64 fGainFMD3i(0),
65 fGainFMD3o(0),
66 fChi2FMD1i(0),
67 fChi2FMD2i(0),
68 fChi2FMD2o(0),
69 fChi2FMD3i(0),
70 fChi2FMD3o(0)
f631f26d 71{
6cf6e7a0 72 // Constructor
73 //
74 // Parameters:
75 // None
427e8f99 76 fCurrentPulse.Reset(0);
77 fCurrentChannel.Reset(0);
f631f26d 78 fOutputFile.open("gains.csv");
a31ea3ce 79 fDiagnosticsFilename = "diagnosticsGain.root";
f631f26d 80}
81
82//_____________________________________________________________________
cf291b42 83AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA)
2a082c96 84 : AliFMDBaseDA(gainDA),
2a082c96 85 fHighPulse(gainDA.fHighPulse),
86 fEventsPerChannel(gainDA.fEventsPerChannel),
87 fCurrentPulse(gainDA.fCurrentPulse),
88 fCurrentChannel(gainDA.fCurrentChannel),
89 fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
90 fSummaryGains(gainDA.fSummaryGains),
91 fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip),
92 fGainFMD1i(gainDA.fGainFMD1i),
93 fGainFMD2i(gainDA.fGainFMD2i),
94 fGainFMD2o(gainDA.fGainFMD2o),
95 fGainFMD3i(gainDA.fGainFMD3i),
96 fGainFMD3o(gainDA.fGainFMD3o),
97 fChi2FMD1i(gainDA.fChi2FMD1i),
98 fChi2FMD2i(gainDA.fChi2FMD2i),
99 fChi2FMD2o(gainDA.fChi2FMD2o),
100 fChi2FMD3i(gainDA.fChi2FMD3i),
101 fChi2FMD3o(gainDA.fChi2FMD3o)
cf291b42 102{
6cf6e7a0 103 // Copy Constructor
104 //
105 // Parameters:
106 // gainDA Object to copy from
427e8f99 107 fCurrentPulse.Reset(0);
108 fCurrentChannel.Reset(0);
f631f26d 109}
110
111//_____________________________________________________________________
cf291b42 112AliFMDGainDA::~AliFMDGainDA()
113{
6cf6e7a0 114 // Destructor
115 //
116 // Parameters:
117 // None
f631f26d 118}
119
120//_____________________________________________________________________
cf291b42 121void AliFMDGainDA::Init()
122{
6cf6e7a0 123 // Initialize
124 //
125 // Parameters:
126 // None
427e8f99 127 Int_t nEventsRequired = 0;
22017a7e 128
6cf6e7a0 129 for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++) {
130 Int_t nEvents = 0;
131 if(fPulseSize.At(idx))
132 nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
133 fEventsPerChannel.AddAt(nEvents,idx);
134 if(nEvents>nEventsRequired)
a31ea3ce 135 nEventsRequired = nEvents * fNumberOfStripsPerChip;
6cf6e7a0 136 }
427e8f99 137 SetRequiredEvents(nEventsRequired);
a31ea3ce 138}
139
140//_____________________________________________________________________
141void AliFMDGainDA::InitContainer(TDirectory* dir)
142{
143 AliFMDBaseDA::InitContainer(dir);
f631f26d 144
145 for(UShort_t det=1;det<=3;det++) {
a31ea3ce 146 UShort_t sr = (det == 1 ? 1 : 0);
147 for (UShort_t ir = sr; ir < 2; ir++) {
f631f26d 148 Char_t ring = (ir == 0 ? 'O' : 'I');
149 UShort_t nsec = (ir == 0 ? 40 : 20);
a31ea3ce 150 UShort_t nva = (ir == 0 ? 2 : 4);
f631f26d 151 for(UShort_t sec =0; sec < nsec; sec++) {
3490bd31 152 Array* sectorArray = GetSectorArray(det, ring, sec);
153 Array* cache = new Array(nva);
a31ea3ce 154 cache->SetName("Cache");
155 cache->SetOwner();
156 Int_t n = sectorArray->GetEntriesFast();
157 sectorArray->AddAtAndExpand(cache, n);
158 for(UShort_t va = 0; va < nva; va++) {
159 TH1S* hChannel = new TH1S(Form("FMD%d%c[%02d]_va%d",
160 det,ring,sec,va),
161 Form("FMD%d%c[%02d] VA%d Cache",
162 det,ring,sec,va),
163 1024,-.5,1023.5);
f631f26d 164 hChannel->SetDirectory(0);
a31ea3ce 165 cache->AddAtAndExpand(hChannel,va);
f631f26d 166 }
167 }
168 }
169 }
170}
171
172//_____________________________________________________________________
3490bd31 173void AliFMDGainDA::AddChannelContainer(Array* stripArray,
cf291b42 174 UShort_t det ,
175 Char_t ring,
f631f26d 176 UShort_t sec,
cf291b42 177 UShort_t strip)
178{
6cf6e7a0 179 // Make a channel container
180 //
181 // Parameters:
182 // sectorArray Sectors
183 // det Detector number
184 // ring Ring identifier
185 // sec Sector number
186 // strip Strip number
a31ea3ce 187
188 AliFMDParameters* pars = AliFMDParameters::Instance();
189 UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
190 Int_t halfring = GetHalfringIndex(det,ring,board/16);
191 Int_t dPulse = fPulseSize.At(halfring);
192 Int_t nPulses = dPulse > 0 ? fHighPulse / dPulse : 0;
193
194 TGraphErrors* gChannel = new TGraphErrors(nPulses);
195 gChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
196 gChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC",
cf291b42 197 det, ring, sec, strip));
a31ea3ce 198 stripArray->AddAtAndExpand(gChannel,0);
199}
200
201//_____________________________________________________________________
3490bd31 202void AliFMDGainDA::AddSectorSummary(Array* sectorArray,
a31ea3ce 203 UShort_t det,
204 Char_t ring,
205 UShort_t sec,
206 UShort_t nStr)
207{
208 TH1F* summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]",
209 det, ring, sec),
210 nStr, -.5, nStr-.5);
211 summary->SetXTitle("Strip");
212 summary->SetYTitle("Gain [ADC/DAC]");
213 summary->SetDirectory(0);
214
215 Int_t n = sectorArray->GetEntriesFast();
216 sectorArray->AddAtAndExpand(summary, n);
f631f26d 217}
218
219//_____________________________________________________________________
6cf6e7a0 220void AliFMDGainDA::FillChannels(AliFMDDigit* digit)
221{
222 // Fill data into histogram
223 //
224 // Parameters:
225 // digit Digit to get the data from
f631f26d 226
227 UShort_t det = digit->Detector();
228 Char_t ring = digit->Ring();
229 UShort_t sec = digit->Sector();
230 UShort_t strip = digit->Strip();
231
95779ff5 232 //Strip is always seen as the first in a VA chip. All other strips are junk.
233 //Strips are counted from zero on even sectors and from 511 on odd sectors...
234
235 if((sec%2) && ((strip+1) % fNumberOfStripsPerChip)) return;
236 if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
237
a31ea3ce 238 UShort_t vaChip = strip / fNumberOfStripsPerChip;
239 TH1S* hChannel = GetChannelHistogram(det, ring, sec, vaChip);
f631f26d 240 hChannel->Fill(digit->Counts());
241 UpdatePulseAndADC(det,ring,sec,strip);
242}
243
244//_____________________________________________________________________
2a082c96 245void AliFMDGainDA::MakeSummary(UShort_t det, Char_t ring)
246{
8cc1cb63 247 //
248 // Create summary hists for FMD gains and chi2 of the fits
249 //
2a082c96 250 switch (det) {
251 case 1:
252 fGainFMD1i = MakeSummaryHistogram("gain", "Gains", det, ring);
253 fChi2FMD1i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
254 break;
255 case 2:
256 switch (ring) {
257 case 'I': case 'i':
258 fGainFMD2i = MakeSummaryHistogram("gain", "Gains", det, ring);
259 fChi2FMD2i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
260 break;
261 case 'O': case 'o':
262 fGainFMD2o = MakeSummaryHistogram("gain", "Gains", det, ring);
263 fChi2FMD2o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
264 break;
265 }
266 break;
267 case 3:
268 switch (ring) {
269 case 'I': case 'i':
270 fGainFMD3i = MakeSummaryHistogram("gain", "Gains", det, ring);
271 fChi2FMD3i = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
272 break;
273 case 'O': case 'o':
274 fGainFMD3o = MakeSummaryHistogram("gain", "Gains", det, ring);
275 fChi2FMD3o = MakeSummaryHistogram("chi2", "#Chi^{2}/NDF", det, ring);
276 break;
277 }
278 break;
279 }
280}
281
282//_____________________________________________________________________
f631f26d 283void AliFMDGainDA::Analyse(UShort_t det,
e9c06036 284 Char_t ring,
f631f26d 285 UShort_t sec,
6cf6e7a0 286 UShort_t strip)
287{
288 // Analyse result of a single strip
289 //
290 // Parameters:
291 // det Detector number
292 // ring Ring identifier
293 // sec Sector number
294 // strip Strip number
f631f26d 295 TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
296 if(!grChannel->GetN()) {
95779ff5 297 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
298 det, ring , sec, strip));
f631f26d 299 return;
300 }
f8022830 301 TF1* fitFunc = new TF1("fitFunc","pol1",-10,280);
302 fitFunc->SetParameters(100,3);
303 grChannel->Fit(fitFunc,"Q0","",0,fHighPulse);
304 if (grChannel->GetListOfFunctions())
305 grChannel->GetListOfFunctions()->Remove(fitFunc);
306 gROOT->GetListOfFunctions()->Remove(fitFunc);
307
f7f0b643 308
4fa150a7 309 Float_t gain = -1;
310 Float_t error = -1;
311 Float_t chi2ndf = -1;
f8022830 312 if((fitFunc->GetParameter(1)) == (fitFunc->GetParameter(1))) {
313 gain = fitFunc->GetParameter(1);
314 error = fitFunc->GetParError(1);
315 if(fitFunc->GetNDF())
316 chi2ndf = fitFunc->GetChisquare() / fitFunc->GetNDF();
4fa150a7 317 }
f8022830 318
f7f0b643 319 fOutputFile << det << ','
320 << ring << ','
321 << sec << ','
322 << strip << ','
4fa150a7 323 << gain << ','
324 << error << ','
f631f26d 325 << chi2ndf <<"\n";
326
22017a7e 327 //due to RCU trouble, first strips on VAs are excluded
276b1261 328 // if(strip%128 != 0) {
22017a7e 329
f8022830 330 fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc->GetParameter(1));
331 fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc->GetParError(1));
276b1261 332
333 fCurrentSummaryStrip++;
2a082c96 334
335 TH2* hGain = 0;
336 TH2* hChi2 = 0;
337 switch (det) {
338 case 1: hGain = fGainFMD1i; hChi2 = fChi2FMD1i; break;
339 case 2:
340 switch (ring) {
341 case 'I': hGain = fGainFMD2i; hChi2 = fChi2FMD2i; break;
342 case 'O': hGain = fGainFMD2o; hChi2 = fChi2FMD2o; break;
343 }
344 break;
345 case 3:
346 switch (ring) {
347 case 'I': hGain = fGainFMD3i; hChi2 = fChi2FMD3i; break;
348 case 'O': hGain = fGainFMD3o; hChi2 = fChi2FMD3o; break;
349 }
350 break;
351 }
352 if (hGain && hChi2) {
353 Int_t bin = hGain->FindBin(sec, strip);
354 hGain->SetBinContent(bin, gain);
355 hGain->SetBinError(bin, error);
356 hChi2->SetBinContent(bin, chi2ndf);
357 }
358
276b1261 359 // }
f631f26d 360 if(fSaveHistograms) {
a31ea3ce 361 TH1F* summary = GetSectorSummary(det, ring, sec);
f8022830 362 summary->SetBinContent(strip+1, fitFunc->GetParameter(1));
363 summary->SetBinError(strip+1, fitFunc->GetParError(1));
e9c06036 364 }
f8022830 365 delete fitFunc;
f631f26d 366}
367
368//_____________________________________________________________________
3bae5d02 369void AliFMDGainDA::Terminate(TFile* diagFile)
370{
6cf6e7a0 371 // End of file
372 //
373 // Parameters:
374 // None
7bce699b 375 if(diagFile) {
376 diagFile->cd();
377 fSummaryGains.Write();
378 }
3bae5d02 379}
380
381//_____________________________________________________________________
e9c06036 382void AliFMDGainDA::WriteHeaderToFile()
383{
6cf6e7a0 384 // Write header to the output file
385 //
386 // Parameters:
387 // None
9c958f2b 388 AliFMDParameters* pars = AliFMDParameters::Instance();
389 fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
6cf6e7a0 390 TDatime now;
391 fOutputFile << "# This file created from run # " << fRunno
392 << " @ " << now.AsString() << std::endl;
f7f0b643 393 fOutputFile.write("# Detector, "
394 "Ring, "
395 "Sector, "
e9c06036 396 "Strip, "
397 "Gain, "
398 "Error, "
a31ea3ce 399 "Chi2/NDF \n",56);
f631f26d 400}
401
402//_____________________________________________________________________
403TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
e9c06036 404 Char_t ring,
f631f26d 405 UShort_t sec,
a31ea3ce 406 UShort_t va)
e9c06036 407{
6cf6e7a0 408 // Get the current histogram of a single strip
409 //
410 // Parameters:
411 // det Detector number
412 // ring Ring identifier
413 // sec Sector number
a31ea3ce 414 // va VA number
3490bd31 415 Array* secArray = GetSectorArray(det, ring, sec);
a31ea3ce 416 Int_t n = secArray->GetEntriesFast();
3490bd31 417 Array* cache = static_cast<Array*>(secArray->At(n-1));
a31ea3ce 418 TH1S* hChannel = static_cast<TH1S*>(cache->At(va));
f631f26d 419
420 return hChannel;
421}
422
423//_____________________________________________________________________
cf291b42 424TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
e9c06036 425 Char_t ring,
cf291b42 426 UShort_t sec,
e9c06036 427 UShort_t strip)
428{
6cf6e7a0 429 // Get the graph of a single strip
430 //
431 // Parameters:
432 // det Detector number
433 // ring Ring identifier
434 // sec Sector number
435 // strip Strip number
3490bd31 436 Array* stripArray = GetStripArray(det, ring, sec, strip);
a31ea3ce 437 TGraphErrors* hChannel = static_cast<TGraphErrors*>(stripArray->At(0));
f631f26d 438
439 return hChannel;
440}
441
442//_____________________________________________________________________
a31ea3ce 443TH1F* AliFMDGainDA::GetSectorSummary(UShort_t det,
444 Char_t ring,
445 UShort_t sec)
446{
3490bd31 447 Array* secArray = GetSectorArray(det, ring, sec);
a31ea3ce 448 Int_t n = secArray->GetEntriesFast();
449 return static_cast<TH1F*>(secArray->At(n-2)); // Cache added later
450}
451
452//_____________________________________________________________________
cf291b42 453void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
454 Char_t ring,
455 UShort_t sec,
456 UShort_t strip)
457{
6cf6e7a0 458 // Go to next pulse size
459 //
460 // Parameters:
461 // det Detector number
462 // ring Ring identifier
463 // sec Sector number
464 // strip Strip number
f7f0b643 465
a31ea3ce 466 AliFMDParameters* pars = AliFMDParameters::Instance();
467 UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
468 Int_t halfring = GetHalfringIndex(det,ring,board/16);
22017a7e 469
f7f0b643 470 if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
aea2699c 471 return;
95779ff5 472
a31ea3ce 473 if ((sec % 2) && ((strip+1) % fNumberOfStripsPerChip)) return;
474 if (((sec + 1) % 2) && (strip % fNumberOfStripsPerChip)) return;
95779ff5 475
b995fc28 476 if(((GetCurrentEvent()) % fPulseLength.At(halfring))
477 && GetCurrentEvent() > 0) return;
427e8f99 478
a31ea3ce 479 Int_t vaChip = strip/fNumberOfStripsPerChip;
cf291b42 480 TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
f631f26d 481
482 if(!hChannel->GetEntries()) {
cf291b42 483 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
484 det, ring , sec, strip));
f631f26d 485 return;
486 }
487 Double_t mean = hChannel->GetMean();
488 Double_t rms = hChannel->GetRMS();
b995fc28 489 Double_t pulse = (Double_t(fCurrentPulse.At(halfring))
490 * fPulseSize.At(halfring));
cf291b42 491 Int_t firstBin = hChannel->GetXaxis()->GetFirst();
492 Int_t lastBin = hChannel->GetXaxis()->GetLast();
9c958f2b 493 hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
f631f26d 494
cf291b42 495 mean = hChannel->GetMean();
496 rms = hChannel->GetRMS();
efd39294 497
498 hChannel->GetXaxis()->SetRange(firstBin,lastBin);
499
a31ea3ce 500 Int_t channelOff = ((GetCurrentEvent()-1)
501 / ((fPulseLength.At(halfring)*fHighPulse)
502 / fPulseSize.At(halfring)));
503 Int_t channelNo = (strip + ((sec % 2 == 1) ? -1 : 1) * channelOff);
504 TGraphErrors* channel = GetChannel(det,ring,sec,channelNo);
505
427e8f99 506 channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
507 channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
f631f26d 508
509 if(fSaveHistograms) {
a31ea3ce 510 TH1S* out =
511 static_cast<TH1S*>(hChannel->Clone(Form("FMD%d%c[%02d,%03d]_0x%02x",
512 det, ring, sec, channelNo,
513 int(pulse))));
514 out->SetTitle(Form("FMD%d%c[%02d,%03d] DAC=0x%02x (%3d)",
515 det, ring, sec, channelNo, int(pulse), int(pulse)));
516 out->SetDirectory(0);
3490bd31 517 Array* arr = GetStripArray(det, ring, sec, channelNo);
a31ea3ce 518 arr->AddAtAndExpand(out, fCurrentPulse.At(halfring)+1);
f631f26d 519 }
22017a7e 520
f631f26d 521 hChannel->Reset();
522
523}
524
525//_____________________________________________________________________
cf291b42 526void AliFMDGainDA::ResetPulseAndUpdateChannel()
527{
6cf6e7a0 528 // Reset all
529 //
530 // Parameters:
531 // None
427e8f99 532 fCurrentPulse.Reset(0);
f631f26d 533}
534
535//_____________________________________________________________________
cf291b42 536void AliFMDGainDA::FinishEvent()
537{
6cf6e7a0 538 // End of event
539 //
540 // Parameters:
541 // None
22017a7e 542 for(UShort_t det=1; det<=3;det++) {
543 UShort_t firstring = (det == 1 ? 1 : 0);
544 for(UShort_t iring = firstring; iring <=1;iring++) {
545 Char_t ring = (iring == 1 ? 'I' : 'O');
546 for(UShort_t board =0 ; board <=1; board++) {
547 Int_t idx = GetHalfringIndex(det,ring,board);
548
6cf6e7a0 549 if(!fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
22017a7e 550 continue;
6cf6e7a0 551 if(GetCurrentEvent()>0 &&
552 ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
22017a7e 553 fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
554
6cf6e7a0 555 if(GetCurrentEvent()>0 &&
556 ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
22017a7e 557 fCurrentPulse.AddAt(0,idx);
558 }
559 }
427e8f99 560 }
f631f26d 561}
562//_____________________________________________________________________
563//
564//EOF
565//