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