]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - FMD/AliFMDGainDA.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Updates
[u/mrichter/AliRoot.git] / FMD / AliFMDGainDA.cxx
... / ...
CommitLineData
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 "AliFMDAltroMapping.h"
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>
41#include <TDatime.h>
42#include <TH2.h>
43
44//_____________________________________________________________________
45ClassImp(AliFMDGainDA)
46#if 0 // Do not delete - here to let Emacs indent properly
47;
48#endif
49
50//_____________________________________________________________________
51AliFMDGainDA::AliFMDGainDA()
52 : AliFMDBaseDA(),
53 fHighPulse(256),
54 fEventsPerChannel(10),
55 fCurrentPulse(10),
56 fCurrentChannel(10),
57 fNumberOfStripsPerChip(128),
58 fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
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)
70{
71 // Constructor
72 //
73 // Parameters:
74 // None
75 fCurrentPulse.Reset(0);
76 fCurrentChannel.Reset(0);
77 fOutputFile.open("gains.csv");
78 fDiagnosticsFilename = "diagnosticsGain.root";
79}
80
81//_____________________________________________________________________
82AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA)
83 : AliFMDBaseDA(gainDA),
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)
101{
102 // Copy Constructor
103 //
104 // Parameters:
105 // gainDA Object to copy from
106 fCurrentPulse.Reset(0);
107 fCurrentChannel.Reset(0);
108}
109
110//_____________________________________________________________________
111AliFMDGainDA::~AliFMDGainDA()
112{
113 // Destructor
114 //
115 // Parameters:
116 // None
117}
118
119//_____________________________________________________________________
120void AliFMDGainDA::Init()
121{
122 // Initialize
123 //
124 // Parameters:
125 // None
126 Int_t nEventsRequired = 0;
127
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)
134 nEventsRequired = nEvents * fNumberOfStripsPerChip;
135 }
136 SetRequiredEvents(nEventsRequired);
137}
138
139//_____________________________________________________________________
140void AliFMDGainDA::InitContainer(TDirectory* dir)
141{
142 AliFMDBaseDA::InitContainer(dir);
143
144 for(UShort_t det=1;det<=3;det++) {
145 UShort_t sr = (det == 1 ? 1 : 0);
146 for (UShort_t ir = sr; ir < 2; ir++) {
147 Char_t ring = (ir == 0 ? 'O' : 'I');
148 UShort_t nsec = (ir == 0 ? 40 : 20);
149 UShort_t nva = (ir == 0 ? 2 : 4);
150 for(UShort_t sec =0; sec < nsec; sec++) {
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);
163 hChannel->SetDirectory(0);
164 cache->AddAtAndExpand(hChannel,va);
165 }
166 }
167 }
168 }
169}
170
171//_____________________________________________________________________
172void AliFMDGainDA::AddChannelContainer(TObjArray* stripArray,
173 UShort_t det ,
174 Char_t ring,
175 UShort_t sec,
176 UShort_t strip)
177{
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
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",
196 det, ring, sec, strip));
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);
216}
217
218//_____________________________________________________________________
219void AliFMDGainDA::FillChannels(AliFMDDigit* digit)
220{
221 // Fill data into histogram
222 //
223 // Parameters:
224 // digit Digit to get the data from
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
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
237 UShort_t vaChip = strip / fNumberOfStripsPerChip;
238 TH1S* hChannel = GetChannelHistogram(det, ring, sec, vaChip);
239 hChannel->Fill(digit->Counts());
240 UpdatePulseAndADC(det,ring,sec,strip);
241}
242
243//_____________________________________________________________________
244void AliFMDGainDA::MakeSummary(UShort_t det, Char_t ring)
245{
246 //
247 // Create summary hists for FMD gains and chi2 of the fits
248 //
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//_____________________________________________________________________
282void AliFMDGainDA::Analyse(UShort_t det,
283 Char_t ring,
284 UShort_t sec,
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
294 TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
295 if(!grChannel->GetN()) {
296 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
297 det, ring , sec, strip));
298 return;
299 }
300 TF1 fitFunc("fitFunc","pol1",-10,280);
301 fitFunc.SetParameters(100,3);
302 grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
303
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
314 fOutputFile << det << ','
315 << ring << ','
316 << sec << ','
317 << strip << ','
318 << gain << ','
319 << error << ','
320 << chi2ndf <<"\n";
321
322 //due to RCU trouble, first strips on VAs are excluded
323 // if(strip%128 != 0) {
324
325 fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
326 fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
327
328 fCurrentSummaryStrip++;
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
354 // }
355 if(fSaveHistograms) {
356 TH1F* summary = GetSectorSummary(det, ring, sec);
357 summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
358 summary->SetBinError(strip+1, fitFunc.GetParError(1));
359 }
360}
361
362//_____________________________________________________________________
363void AliFMDGainDA::Terminate(TFile* diagFile)
364{
365 // End of file
366 //
367 // Parameters:
368 // None
369 if(diagFile) {
370 diagFile->cd();
371 fSummaryGains.Write();
372 }
373}
374
375//_____________________________________________________________________
376void AliFMDGainDA::WriteHeaderToFile()
377{
378 // Write header to the output file
379 //
380 // Parameters:
381 // None
382 AliFMDParameters* pars = AliFMDParameters::Instance();
383 fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
384 TDatime now;
385 fOutputFile << "# This file created from run # " << fRunno
386 << " @ " << now.AsString() << std::endl;
387 fOutputFile.write("# Detector, "
388 "Ring, "
389 "Sector, "
390 "Strip, "
391 "Gain, "
392 "Error, "
393 "Chi2/NDF \n",56);
394}
395
396//_____________________________________________________________________
397TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
398 Char_t ring,
399 UShort_t sec,
400 UShort_t va)
401{
402 // Get the current histogram of a single strip
403 //
404 // Parameters:
405 // det Detector number
406 // ring Ring identifier
407 // sec Sector number
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));
413
414 return hChannel;
415}
416
417//_____________________________________________________________________
418TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
419 Char_t ring,
420 UShort_t sec,
421 UShort_t strip)
422{
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
430 TObjArray* stripArray = GetStripArray(det, ring, sec, strip);
431 TGraphErrors* hChannel = static_cast<TGraphErrors*>(stripArray->At(0));
432
433 return hChannel;
434}
435
436//_____________________________________________________________________
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//_____________________________________________________________________
447void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
448 Char_t ring,
449 UShort_t sec,
450 UShort_t strip)
451{
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
459
460 AliFMDParameters* pars = AliFMDParameters::Instance();
461 UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
462 Int_t halfring = GetHalfringIndex(det,ring,board/16);
463
464 if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
465 return;
466
467 if ((sec % 2) && ((strip+1) % fNumberOfStripsPerChip)) return;
468 if (((sec + 1) % 2) && (strip % fNumberOfStripsPerChip)) return;
469
470 if(((GetCurrentEvent()) % fPulseLength.At(halfring))
471 && GetCurrentEvent() > 0) return;
472
473 Int_t vaChip = strip/fNumberOfStripsPerChip;
474 TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
475
476 if(!hChannel->GetEntries()) {
477 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
478 det, ring , sec, strip));
479 return;
480 }
481 Double_t mean = hChannel->GetMean();
482 Double_t rms = hChannel->GetRMS();
483 Double_t pulse = (Double_t(fCurrentPulse.At(halfring))
484 * fPulseSize.At(halfring));
485 Int_t firstBin = hChannel->GetXaxis()->GetFirst();
486 Int_t lastBin = hChannel->GetXaxis()->GetLast();
487 hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
488
489 mean = hChannel->GetMean();
490 rms = hChannel->GetRMS();
491
492 hChannel->GetXaxis()->SetRange(firstBin,lastBin);
493
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
500 channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
501 channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
502
503 if(fSaveHistograms) {
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);
513 }
514
515 hChannel->Reset();
516
517}
518
519//_____________________________________________________________________
520void AliFMDGainDA::ResetPulseAndUpdateChannel()
521{
522 // Reset all
523 //
524 // Parameters:
525 // None
526 fCurrentPulse.Reset(0);
527}
528
529//_____________________________________________________________________
530void AliFMDGainDA::FinishEvent()
531{
532 // End of event
533 //
534 // Parameters:
535 // None
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
543 if(!fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
544 continue;
545 if(GetCurrentEvent()>0 &&
546 ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
547 fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
548
549 if(GetCurrentEvent()>0 &&
550 ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
551 fCurrentPulse.AddAt(0,idx);
552 }
553 }
554 }
555}
556//_____________________________________________________________________
557//
558//EOF
559//