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