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