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