]> 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
Code to draw calibrations
[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(),
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{
8cc1cb63 223 //
224 // Create summary hists for FMD gains and chi2 of the fits
225 //
2a082c96 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
f631f26d 258//_____________________________________________________________________
259void AliFMDGainDA::Analyse(UShort_t det,
e9c06036 260 Char_t ring,
f631f26d 261 UShort_t sec,
6cf6e7a0 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
f631f26d 271 TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
272 if(!grChannel->GetN()) {
95779ff5 273 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
274 det, ring , sec, strip));
f631f26d 275 return;
276 }
9c958f2b 277 TF1 fitFunc("fitFunc","pol1",-10,280);
278 fitFunc.SetParameters(100,3);
e9c06036 279 grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
f7f0b643 280
4fa150a7 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
f7f0b643 291 fOutputFile << det << ','
292 << ring << ','
293 << sec << ','
294 << strip << ','
4fa150a7 295 << gain << ','
296 << error << ','
f631f26d 297 << chi2ndf <<"\n";
298
22017a7e 299 //due to RCU trouble, first strips on VAs are excluded
276b1261 300 // if(strip%128 != 0) {
22017a7e 301
276b1261 302 fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
303 fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
304
305 fCurrentSummaryStrip++;
2a082c96 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
276b1261 331 // }
f631f26d 332 if(fSaveHistograms) {
e9c06036 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 }
f631f26d 354}
355
356//_____________________________________________________________________
3bae5d02 357void AliFMDGainDA::Terminate(TFile* diagFile)
358{
6cf6e7a0 359 // End of file
360 //
361 // Parameters:
362 // None
7bce699b 363 if(diagFile) {
364 diagFile->cd();
365 fSummaryGains.Write();
366 }
3bae5d02 367}
368
369//_____________________________________________________________________
e9c06036 370void AliFMDGainDA::WriteHeaderToFile()
371{
6cf6e7a0 372 // Write header to the output file
373 //
374 // Parameters:
375 // None
9c958f2b 376 AliFMDParameters* pars = AliFMDParameters::Instance();
377 fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
6cf6e7a0 378 TDatime now;
379 fOutputFile << "# This file created from run # " << fRunno
380 << " @ " << now.AsString() << std::endl;
f7f0b643 381 fOutputFile.write("# Detector, "
382 "Ring, "
383 "Sector, "
e9c06036 384 "Strip, "
385 "Gain, "
386 "Error, "
f7f0b643 387 "Chi2/NDF \n",56);
f631f26d 388
389}
390
391//_____________________________________________________________________
392TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
e9c06036 393 Char_t ring,
f631f26d 394 UShort_t sec,
e9c06036 395 UShort_t strip)
396{
6cf6e7a0 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
f631f26d 404
6cf6e7a0 405 UShort_t lRing = 1;
f631f26d 406 if(ring == 'O')
6cf6e7a0 407 lRing = 0;
f631f26d 408
409
410 TObjArray* detArray = static_cast<TObjArray*>(fGainArray.At(det));
6cf6e7a0 411 TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(lRing));
f631f26d 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//_____________________________________________________________________
cf291b42 419TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
e9c06036 420 Char_t ring,
cf291b42 421 UShort_t sec,
e9c06036 422 UShort_t strip)
423{
6cf6e7a0 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
cf291b42 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));
f631f26d 436
437 return hChannel;
438}
439
440//_____________________________________________________________________
cf291b42 441void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
442 Char_t ring,
443 UShort_t sec,
444 UShort_t strip)
445{
6cf6e7a0 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
f7f0b643 453
427e8f99 454 AliFMDParameters* pars = AliFMDParameters::Instance();
b995fc28 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(
22017a7e 459 Int_t halfring = GetHalfringIndex(det,ring,board/16);
460
f7f0b643 461 if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
aea2699c 462 return;
95779ff5 463
464 if((sec%2) && ((strip+1) % fNumberOfStripsPerChip)) return;
465
466 if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
467
b995fc28 468 if(((GetCurrentEvent()) % fPulseLength.At(halfring))
469 && GetCurrentEvent() > 0) return;
427e8f99 470
cf291b42 471 Int_t vaChip = strip/fNumberOfStripsPerChip;
472 TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
f631f26d 473
474 if(!hChannel->GetEntries()) {
cf291b42 475 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
476 det, ring , sec, strip));
f631f26d 477 return;
478 }
479 Double_t mean = hChannel->GetMean();
480 Double_t rms = hChannel->GetRMS();
b995fc28 481 Double_t pulse = (Double_t(fCurrentPulse.At(halfring))
482 * fPulseSize.At(halfring));
cf291b42 483 Int_t firstBin = hChannel->GetXaxis()->GetFirst();
484 Int_t lastBin = hChannel->GetXaxis()->GetLast();
9c958f2b 485 hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
f631f26d 486
cf291b42 487 mean = hChannel->GetMean();
488 rms = hChannel->GetRMS();
efd39294 489
490 hChannel->GetXaxis()->SetRange(firstBin,lastBin);
491
95779ff5 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)));
f631f26d 501
502 TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
503
427e8f99 504 channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
505 channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
f631f26d 506
507 if(fSaveHistograms) {
f7f0b643 508 gDirectory->cd(GetStripPath(det,ring,sec,channelNumber));
509 hChannel->Write(Form("%s_pulse_%03d",hChannel->GetName(),(Int_t)pulse));
510
f631f26d 511 }
22017a7e 512
f631f26d 513 hChannel->Reset();
514
515}
516
517//_____________________________________________________________________
cf291b42 518void AliFMDGainDA::ResetPulseAndUpdateChannel()
519{
6cf6e7a0 520 // Reset all
521 //
522 // Parameters:
523 // None
427e8f99 524 fCurrentPulse.Reset(0);
f631f26d 525}
526
527//_____________________________________________________________________
cf291b42 528void AliFMDGainDA::FinishEvent()
529{
6cf6e7a0 530 // End of event
531 //
532 // Parameters:
533 // None
22017a7e 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
6cf6e7a0 541 if(!fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
22017a7e 542 continue;
6cf6e7a0 543 if(GetCurrentEvent()>0 &&
544 ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
22017a7e 545 fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
546
6cf6e7a0 547 if(GetCurrentEvent()>0 &&
548 ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
22017a7e 549 fCurrentPulse.AddAt(0,idx);
550 }
551 }
427e8f99 552 }
f631f26d 553}
554//_____________________________________________________________________
555//
556//EOF
557//