Connection to tree in Notify() seems to be safer.
[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"
32#include "iostream"
33#include "fstream"
34#include "AliLog.h"
35#include "TF1.h"
36#include "TH1.h"
37#include "TMath.h"
38#include "TGraphErrors.h"
39#include "AliFMDParameters.h"
b995fc28 40#include "AliFMDAltroMapping.h"
f631f26d 41
42//_____________________________________________________________________
43ClassImp(AliFMDGainDA)
cf291b42 44#if 0 // Do not delete - here to let Emacs indent properly
45;
46#endif
f631f26d 47
48//_____________________________________________________________________
cf291b42 49AliFMDGainDA::AliFMDGainDA()
50 : AliFMDBaseDA(),
51 fGainArray(),
cf291b42 52 fHighPulse(256),
22017a7e 53 fEventsPerChannel(10),
54 fCurrentPulse(10),
55 fCurrentChannel(10),
3bae5d02 56 fNumberOfStripsPerChip(128),
57 fSummaryGains("GainsSummary","Summary of gains",51200,0,51200),
f14ede67 58 fCurrentSummaryStrip(1)
f631f26d 59{
427e8f99 60 fCurrentPulse.Reset(0);
61 fCurrentChannel.Reset(0);
f631f26d 62 fOutputFile.open("gains.csv");
cf291b42 63 fGainArray.SetOwner();
f631f26d 64}
65
66//_____________________________________________________________________
cf291b42 67AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA)
68 : AliFMDBaseDA(gainDA),
69 fGainArray(gainDA.fGainArray),
cf291b42 70 fHighPulse(gainDA.fHighPulse),
cf291b42 71 fEventsPerChannel(gainDA.fEventsPerChannel),
72 fCurrentPulse(gainDA.fCurrentPulse),
73 fCurrentChannel(gainDA.fCurrentChannel),
3bae5d02 74 fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip),
75 fSummaryGains(gainDA.fSummaryGains),
76 fCurrentSummaryStrip(gainDA.fCurrentSummaryStrip)
cf291b42 77{
427e8f99 78 fCurrentPulse.Reset(0);
79 fCurrentChannel.Reset(0);
f631f26d 80}
81
82//_____________________________________________________________________
cf291b42 83AliFMDGainDA::~AliFMDGainDA()
84{
f631f26d 85}
86
87//_____________________________________________________________________
cf291b42 88void AliFMDGainDA::Init()
89{
f631f26d 90
427e8f99 91
92 Int_t nEventsRequired = 0;
22017a7e 93
94 //for(UShort_t det=1; det<=3;det++) {
95 // UShort_t firstring = (det == 1 ? 1 : 0);
96 // for(UShort_t iring = firstring; iring <=1;iring++) {
97 // Char_t ring = (iring == 1 ? 'I' : 'O');
98 // for(UShort_t board =0 ; board <=1; board++) {
99 // Int_t idx = GetHalfringIndex(det,ring,board);
100 for(Int_t idx = 0;idx<fEventsPerChannel.GetSize();idx++)
427e8f99 101 {
22017a7e 102
103 Int_t nEvents = 0;
104 if(fPulseSize.At(idx))
105 nEvents = (fPulseLength.At(idx)*fHighPulse) / fPulseSize.At(idx);
106 fEventsPerChannel.AddAt(nEvents,idx);
427e8f99 107 if(nEvents>nEventsRequired) nEventsRequired = nEvents * fNumberOfStripsPerChip;
108
109 }
22017a7e 110 //}
111 // }
54b167db 112
cf291b42 113 //8 pulser values * 128 strips * 100 samples
427e8f99 114
115
116 SetRequiredEvents(nEventsRequired);
22017a7e 117
f631f26d 118 TObjArray* detArray;
119 TObjArray* ringArray;
120 TObjArray* sectorArray;
121
122 for(UShort_t det=1;det<=3;det++) {
123 detArray = new TObjArray();
124 detArray->SetOwner();
125 fGainArray.AddAtAndExpand(detArray,det);
126 for (UShort_t ir = 0; ir < 2; ir++) {
127 Char_t ring = (ir == 0 ? 'O' : 'I');
128 UShort_t nsec = (ir == 0 ? 40 : 20);
129 UShort_t nstr = (ir == 0 ? 2 : 4);
130 ringArray = new TObjArray();
131 ringArray->SetOwner();
132 detArray->AddAtAndExpand(ringArray,ir);
133 for(UShort_t sec =0; sec < nsec; sec++) {
134 sectorArray = new TObjArray();
135 sectorArray->SetOwner();
136 ringArray->AddAtAndExpand(sectorArray,sec);
137 for(UShort_t strip = 0; strip < nstr; strip++) {
cf291b42 138 TH1S* hChannel = new TH1S(Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
139 Form("hFMD%d%c_%d_%d",det,ring,sec,strip),
140 1024,0,1023);
f631f26d 141 hChannel->SetDirectory(0);
142 sectorArray->AddAtAndExpand(hChannel,strip);
143 }
144 }
145 }
146 }
147}
148
149//_____________________________________________________________________
150void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray,
cf291b42 151 UShort_t det ,
152 Char_t ring,
f631f26d 153 UShort_t sec,
cf291b42 154 UShort_t strip)
155{
f631f26d 156 TGraphErrors* hChannel = new TGraphErrors();
cf291b42 157 hChannel->SetName(Form("FMD%d%c[%02d,%03d]", det, ring, sec, strip));
158 hChannel->SetTitle(Form("FMD%d%c[%02d,%03d] ADC vs DAC",
159 det, ring, sec, strip));
f631f26d 160 sectorArray->AddAtAndExpand(hChannel,strip);
161}
162
163//_____________________________________________________________________
164void AliFMDGainDA::FillChannels(AliFMDDigit* digit) {
165
166 UShort_t det = digit->Detector();
167 Char_t ring = digit->Ring();
168 UShort_t sec = digit->Sector();
169 UShort_t strip = digit->Strip();
170
95779ff5 171 //Strip is always seen as the first in a VA chip. All other strips are junk.
172 //Strips are counted from zero on even sectors and from 511 on odd sectors...
173
174 if((sec%2) && ((strip+1) % fNumberOfStripsPerChip)) return;
175 if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
176
cf291b42 177 Int_t vaChip = strip / fNumberOfStripsPerChip;
178 TH1S* hChannel = GetChannelHistogram(det, ring, sec, vaChip);
f631f26d 179 hChannel->Fill(digit->Counts());
180 UpdatePulseAndADC(det,ring,sec,strip);
181}
182
183//_____________________________________________________________________
184void AliFMDGainDA::Analyse(UShort_t det,
e9c06036 185 Char_t ring,
f631f26d 186 UShort_t sec,
187 UShort_t strip) {
188 TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
189 if(!grChannel->GetN()) {
95779ff5 190 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
191 det, ring , sec, strip));
f631f26d 192 return;
193 }
9c958f2b 194 TF1 fitFunc("fitFunc","pol1",-10,280);
195 fitFunc.SetParameters(100,3);
e9c06036 196 grChannel->Fit("fitFunc","Q0+","",0,fHighPulse);
f7f0b643 197
4fa150a7 198 Float_t gain = -1;
199 Float_t error = -1;
200 Float_t chi2ndf = -1;
201 if((fitFunc.GetParameter(1)) == (fitFunc.GetParameter(1))) {
202 gain = fitFunc.GetParameter(1);
203 error = fitFunc.GetParError(1);
204 if(fitFunc.GetNDF())
205 chi2ndf = fitFunc.GetChisquare() / fitFunc.GetNDF();
206 }
207
f7f0b643 208 fOutputFile << det << ','
209 << ring << ','
210 << sec << ','
211 << strip << ','
4fa150a7 212 << gain << ','
213 << error << ','
f631f26d 214 << chi2ndf <<"\n";
215
22017a7e 216 //due to RCU trouble, first strips on VAs are excluded
276b1261 217 // if(strip%128 != 0) {
22017a7e 218
276b1261 219 fSummaryGains.SetBinContent(fCurrentSummaryStrip,fitFunc.GetParameter(1));
220 fSummaryGains.SetBinError(fCurrentSummaryStrip,fitFunc.GetParError(1));
221
222 fCurrentSummaryStrip++;
223 // }
f631f26d 224 if(fSaveHistograms) {
e9c06036 225 gDirectory->cd(GetSectorPath(det,ring, sec, kTRUE));
226
227 TH1F* summary = dynamic_cast<TH1F*>(gDirectory->Get("Summary"));
228 if (!summary) {
229 Int_t nStr = (ring == 'I' ? 512 : 256);
230 summary = new TH1F("Summary", Form("Summary of gains in FMD%d%c[%02d]",
231 det, ring, sec),
232 nStr, -.5, nStr-.5);
233 summary->SetXTitle("Strip");
234 summary->SetYTitle("Gain [ADC/DAC]");
235 summary->SetDirectory(gDirectory);
236 }
237 summary->SetBinContent(strip+1, fitFunc.GetParameter(1));
238 summary->SetBinError(strip+1, fitFunc.GetParError(1));
239
240 gDirectory->cd(GetStripPath(det,ring,sec,strip, kTRUE));
241 grChannel->SetName(Form("FMD%d%c[%02d,%03d]",det,ring,sec,strip));
242 // grChannel->SetDirectory(gDirectory);
243 grChannel->Write();
244 // grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
245 }
f631f26d 246}
247
248//_____________________________________________________________________
3bae5d02 249void AliFMDGainDA::Terminate(TFile* diagFile)
250{
7bce699b 251 if(diagFile) {
252 diagFile->cd();
253 fSummaryGains.Write();
254 }
3bae5d02 255}
256
257//_____________________________________________________________________
e9c06036 258void AliFMDGainDA::WriteHeaderToFile()
259{
9c958f2b 260 AliFMDParameters* pars = AliFMDParameters::Instance();
261 fOutputFile.write(Form("# %s \n",pars->GetGainShuttleID()),9);
f7f0b643 262 fOutputFile.write("# Detector, "
263 "Ring, "
264 "Sector, "
e9c06036 265 "Strip, "
266 "Gain, "
267 "Error, "
f7f0b643 268 "Chi2/NDF \n",56);
f631f26d 269
270}
271
272//_____________________________________________________________________
273TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
e9c06036 274 Char_t ring,
f631f26d 275 UShort_t sec,
e9c06036 276 UShort_t strip)
277{
f631f26d 278
279 UShort_t Ring = 1;
280 if(ring == 'O')
281 Ring = 0;
282
283
284 TObjArray* detArray = static_cast<TObjArray*>(fGainArray.At(det));
285 TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
286 TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
287 TH1S* hChannel = static_cast<TH1S*>(secArray->At(strip));
288
289 return hChannel;
290}
291
292//_____________________________________________________________________
cf291b42 293TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det,
e9c06036 294 Char_t ring,
cf291b42 295 UShort_t sec,
e9c06036 296 UShort_t strip)
297{
cf291b42 298 UShort_t iring = (ring == 'O' ? 0 : 1);
299 TObjArray* detArray = static_cast<TObjArray*>(fDetectorArray.At(det));
300 TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(iring));
301 TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
302 TGraphErrors* hChannel = static_cast<TGraphErrors*>(secArray->At(strip));
f631f26d 303
304 return hChannel;
305}
306
307//_____________________________________________________________________
cf291b42 308void AliFMDGainDA::UpdatePulseAndADC(UShort_t det,
309 Char_t ring,
310 UShort_t sec,
311 UShort_t strip)
312{
f7f0b643 313
427e8f99 314 AliFMDParameters* pars = AliFMDParameters::Instance();
b995fc28 315 // UInt_t ddl, board,chip,ch;
316 UShort_t board = pars->GetAltroMap()->Sector2Board(ring, sec);
317 // pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,ch);
318 /// pars->GetAltroMap()->Strip2Channel(
22017a7e 319 Int_t halfring = GetHalfringIndex(det,ring,board/16);
320
f7f0b643 321 if(GetCurrentEvent()> (fNumberOfStripsPerChip*fEventsPerChannel.At(halfring)))
aea2699c 322 return;
95779ff5 323
324 if((sec%2) && ((strip+1) % fNumberOfStripsPerChip)) return;
325
326 if(((sec+1)%2) && (strip % fNumberOfStripsPerChip)) return;
327
b995fc28 328 if(((GetCurrentEvent()) % fPulseLength.At(halfring))
329 && GetCurrentEvent() > 0) return;
427e8f99 330
cf291b42 331 Int_t vaChip = strip/fNumberOfStripsPerChip;
332 TH1S* hChannel = GetChannelHistogram(det,ring,sec,vaChip);
f631f26d 333
334 if(!hChannel->GetEntries()) {
cf291b42 335 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",
336 det, ring , sec, strip));
f631f26d 337 return;
338 }
339 Double_t mean = hChannel->GetMean();
340 Double_t rms = hChannel->GetRMS();
b995fc28 341 Double_t pulse = (Double_t(fCurrentPulse.At(halfring))
342 * fPulseSize.At(halfring));
cf291b42 343 Int_t firstBin = hChannel->GetXaxis()->GetFirst();
344 Int_t lastBin = hChannel->GetXaxis()->GetLast();
9c958f2b 345 hChannel->GetXaxis()->SetRangeUser(mean-4*rms,mean+4*rms);
f631f26d 346
cf291b42 347 mean = hChannel->GetMean();
348 rms = hChannel->GetRMS();
efd39294 349
350 hChannel->GetXaxis()->SetRange(firstBin,lastBin);
351
95779ff5 352 Int_t channelNumber = (strip +
353 (GetCurrentEvent()-1)
354 / ((fPulseLength.At(halfring)*fHighPulse)
355 / fPulseSize.At(halfring)));
356 if(sec%2)
357 channelNumber = (strip -
358 (GetCurrentEvent()-1)
359 / ((fPulseLength.At(halfring)*fHighPulse)
360 / fPulseSize.At(halfring)));
f631f26d 361
362 TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
363
427e8f99 364 channel->SetPoint(fCurrentPulse.At(halfring),pulse,mean);
365 channel->SetPointError(fCurrentPulse.At(halfring),0,rms);
f631f26d 366
367 if(fSaveHistograms) {
f7f0b643 368 gDirectory->cd(GetStripPath(det,ring,sec,channelNumber));
369 hChannel->Write(Form("%s_pulse_%03d",hChannel->GetName(),(Int_t)pulse));
370
f631f26d 371 }
22017a7e 372
f631f26d 373 hChannel->Reset();
374
375}
376
377//_____________________________________________________________________
cf291b42 378void AliFMDGainDA::ResetPulseAndUpdateChannel()
379{
427e8f99 380 fCurrentPulse.Reset(0);
f631f26d 381}
382
383//_____________________________________________________________________
cf291b42 384void AliFMDGainDA::FinishEvent()
385{
22017a7e 386 for(UShort_t det=1; det<=3;det++) {
387 UShort_t firstring = (det == 1 ? 1 : 0);
388 for(UShort_t iring = firstring; iring <=1;iring++) {
389 Char_t ring = (iring == 1 ? 'I' : 'O');
390 for(UShort_t board =0 ; board <=1; board++) {
391 Int_t idx = GetHalfringIndex(det,ring,board);
392
393 if( !fPulseLength.At(idx) || !fEventsPerChannel.At(idx))
394 continue;
395 if(GetCurrentEvent()>0 && ((GetCurrentEvent() % fPulseLength.At(idx)) == 0))
396 fCurrentPulse.AddAt(fCurrentPulse.At(idx)+1,idx);
397
398 if(GetCurrentEvent()>0 && ((GetCurrentEvent()) % fEventsPerChannel.At(idx)) == 0)
399 fCurrentPulse.AddAt(0,idx);
400 }
401 }
427e8f99 402 }
f631f26d 403}
404//_____________________________________________________________________
405//
406//EOF
407//