This is the gain detector algorithm that will calibrate the FMD.
[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*/
21// This class contains the implementation of the gain detector algorithms (DA) for the FMD.
22// The data is collected in histograms that are reset for each pulse length after the
23// mean and standard deviation are put into a TGraphErrors object. After a certain number of pulses
24// (usually 8) the graph is fitted to a straight line. The gain is then slope of this line as
25// it combines the known pulse and the response of the detector.
26
27
28#include "AliFMDGainDA.h"
29#include "iostream"
30#include "fstream"
31#include "AliLog.h"
32#include "TF1.h"
33#include "TH1.h"
34#include "TMath.h"
35#include "TGraphErrors.h"
36#include "AliFMDParameters.h"
37
38//_____________________________________________________________________
39ClassImp(AliFMDGainDA)
40
41//_____________________________________________________________________
42AliFMDGainDA::AliFMDGainDA() : AliFMDBaseDA(),
43 fGainArray(),
44 fPulseSize(32),
45 fHighPulse(256),
46 fPulseLength(100),
47 fEventsPerChannel(0),
48 fCurrentPulse(0),
49 fCurrentChannel(0),
50 fNumberOfStripsPerChip(128)
51{
52 fOutputFile.open("gains.csv");
53 fGainArray.SetOwner();
54
55}
56
57//_____________________________________________________________________
58AliFMDGainDA::AliFMDGainDA(const AliFMDGainDA & gainDA) :
59 AliFMDBaseDA(gainDA),
60 fGainArray(gainDA.fGainArray),
61 fPulseSize(gainDA.fPulseSize),
62 fHighPulse(gainDA.fHighPulse),
63 fPulseLength(gainDA.fPulseLength),
64 fEventsPerChannel(gainDA.fEventsPerChannel),
65 fCurrentPulse(gainDA.fCurrentPulse),
66 fCurrentChannel(gainDA.fCurrentChannel),
67 fNumberOfStripsPerChip(gainDA.fNumberOfStripsPerChip)
68{
69
70}
71
72//_____________________________________________________________________
73AliFMDGainDA::~AliFMDGainDA() {
74
75
76
77}
78
79//_____________________________________________________________________
80void AliFMDGainDA::Init() {
81 std::cout<<"FMD Gain DA Init"<<std::endl;
82
83 fEventsPerChannel = (fPulseLength*fHighPulse) / fPulseSize ;
84
85 SetRequiredEvents(fEventsPerChannel*fNumberOfStripsPerChip); //8 pulser values * 128 strips * 100 samples
86 TObjArray* detArray;
87 TObjArray* ringArray;
88 TObjArray* sectorArray;
89
90 for(UShort_t det=1;det<=3;det++) {
91 detArray = new TObjArray();
92 detArray->SetOwner();
93 fGainArray.AddAtAndExpand(detArray,det);
94 for (UShort_t ir = 0; ir < 2; ir++) {
95 Char_t ring = (ir == 0 ? 'O' : 'I');
96 UShort_t nsec = (ir == 0 ? 40 : 20);
97 UShort_t nstr = (ir == 0 ? 2 : 4);
98 ringArray = new TObjArray();
99 ringArray->SetOwner();
100 detArray->AddAtAndExpand(ringArray,ir);
101 for(UShort_t sec =0; sec < nsec; sec++) {
102 sectorArray = new TObjArray();
103 sectorArray->SetOwner();
104 ringArray->AddAtAndExpand(sectorArray,sec);
105 for(UShort_t strip = 0; strip < nstr; strip++) {
106 TH1S* hChannel = new TH1S(Form("hFMD%d%c_%d_%d",det,ring,sec,strip),Form("hFMD%d%c_%d_%d",det,ring,sec,strip),1024,0,1023);
107 hChannel->SetDirectory(0);
108 sectorArray->AddAtAndExpand(hChannel,strip);
109 }
110 }
111 }
112 }
113}
114
115//_____________________________________________________________________
116void AliFMDGainDA::AddChannelContainer(TObjArray* sectorArray,
117 UShort_t det,
118 Char_t ring,
119 UShort_t sec,
120 UShort_t strip) {
121
122 TGraphErrors* hChannel = new TGraphErrors();
123 sectorArray->AddAtAndExpand(hChannel,strip);
124}
125
126//_____________________________________________________________________
127void AliFMDGainDA::FillChannels(AliFMDDigit* digit) {
128
129 UShort_t det = digit->Detector();
130 Char_t ring = digit->Ring();
131 UShort_t sec = digit->Sector();
132 UShort_t strip = digit->Strip();
133
134
135 if(strip%fNumberOfStripsPerChip) return;
136 Int_t VAchip = strip / fNumberOfStripsPerChip;
137 TH1S* hChannel = GetChannelHistogram(det,ring,sec,VAchip);
138 hChannel->Fill(digit->Counts());
139 UpdatePulseAndADC(det,ring,sec,strip);
140}
141
142//_____________________________________________________________________
143void AliFMDGainDA::Analyse(UShort_t det,
144 Char_t ring,
145 UShort_t sec,
146 UShort_t strip) {
147 TGraphErrors* grChannel = GetChannel(det,ring,sec,strip);
148 if(!grChannel->GetN()) {
149 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",det, ring , sec, strip));
150 return;
151 }
152 TF1* fitFunc = new TF1("fitFunc","pol1",-10,280);
153 fitFunc->SetParameters(100,3);
154 grChannel->Fit("fitFunc","Q","Q",0,fHighPulse);
155 AliFMDParameters* pars = AliFMDParameters::Instance();
156 UInt_t ddl, board,chip,channel;
157 pars->Detector2Hardware(det,ring,sec,strip,ddl,board,chip,channel);
158
159 Float_t chi2ndf = 0;
160 if(fitFunc->GetNDF())
161 chi2ndf = fitFunc->GetChisquare() / fitFunc->GetNDF();
162 ddl = ddl + kBaseDDL;
163
164 Int_t relStrip = strip%fNumberOfStripsPerChip;
165
166 fOutputFile << ddl << ','
167 << board << ','
168 << chip << ','
169 << channel << ','
170 << relStrip << ','
171 << fitFunc->GetParameter(1) << ','
172 << fitFunc->GetParError(1) << ','
173 << chi2ndf <<"\n";
174
175
176 if(fSaveHistograms) {
177
178 gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,strip));
179 grChannel->Write(Form("grFMD%d%c_%d_%d",det,ring,sec,strip));
180 }
181
182 delete fitFunc;
183 delete grChannel;
184}
185
186//_____________________________________________________________________
187void AliFMDGainDA::WriteHeaderToFile() {
188 fOutputFile.write("# Gains \n",9);
189 fOutputFile.write("# Rcu, Board, Chip, Channel, Strip, Gain, Error, Chi2/NDF \n",59);
190
191}
192
193//_____________________________________________________________________
194TH1S* AliFMDGainDA::GetChannelHistogram(UShort_t det,
195 Char_t ring,
196 UShort_t sec,
197 UShort_t strip) {
198
199 UShort_t Ring = 1;
200 if(ring == 'O')
201 Ring = 0;
202
203
204 TObjArray* detArray = static_cast<TObjArray*>(fGainArray.At(det));
205 TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
206 TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
207 TH1S* hChannel = static_cast<TH1S*>(secArray->At(strip));
208
209 return hChannel;
210}
211
212//_____________________________________________________________________
213TGraphErrors* AliFMDGainDA::GetChannel(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip) {
214
215 UShort_t Ring = 1;
216 if(ring == 'O')
217 Ring = 0;
218
219
220 TObjArray* detArray = static_cast<TObjArray*>(fDetectorArray.At(det));
221 TObjArray* ringArray = static_cast<TObjArray*>(detArray->At(Ring));
222 TObjArray* secArray = static_cast<TObjArray*>(ringArray->At(sec));
223 TGraphErrors* hChannel = static_cast<TGraphErrors*>(secArray->At(strip));
224
225 return hChannel;
226}
227
228//_____________________________________________________________________
229void AliFMDGainDA::UpdatePulseAndADC(UShort_t det, Char_t ring, UShort_t sec, UShort_t strip) {
230
231 if(strip%fNumberOfStripsPerChip) return;
232 if(((GetCurrentEvent())%fPulseLength) && GetCurrentEvent()>0) return;
233
234 Int_t VAchip = strip/fNumberOfStripsPerChip;
235 TH1S* hChannel = GetChannelHistogram(det,ring,sec,VAchip);
236
237 if(!hChannel->GetEntries()) {
238 AliWarning(Form("No entries for FMD%d%c, sector %d, strip %d",det, ring , sec, strip));
239 return;
240 }
241 Double_t mean = hChannel->GetMean();
242 Double_t rms = hChannel->GetRMS();
243 Double_t pulse = (Double_t)fCurrentPulse*fPulseSize;
244
245 for(Int_t i=1;i<hChannel->GetNbinsX();i++)
246 if(hChannel->GetBinCenter(i)>mean+5*rms || hChannel->GetBinCenter(i)<mean-5*rms)
247 hChannel->SetBinContent(i,0);
248
249 mean = hChannel->GetMean();
250 rms = hChannel->GetRMS();
251
252 Int_t channelNumber = strip + (GetCurrentEvent()-1)/800;
253
254 TGraphErrors* channel = GetChannel(det,ring,sec,channelNumber);
255
256 channel->SetPoint(fCurrentPulse,pulse,mean);
257 channel->SetPointError(fCurrentPulse,0,rms);
258
259 if(fSaveHistograms) {
260 gDirectory->cd(Form("%s:FMD%d%c/sector_%d/strip_%d",fDiagnosticsFilename,det,ring,sec,channelNumber));
261 hChannel->Write(Form("hFMD%d%c_%d_%d_pulse_%d",det,ring,sec,channelNumber,fCurrentPulse));
262 }
263
264 hChannel->Reset();
265
266}
267
268//_____________________________________________________________________
269void AliFMDGainDA::ResetPulseAndUpdateChannel() {
270
271 fCurrentPulse = 0;
272
273}
274
275//_____________________________________________________________________
276void AliFMDGainDA::FinishEvent() {
277
278 if(GetCurrentEvent()>0 && (GetCurrentEvent()%fPulseLength==0))
279 fCurrentPulse++;
280
281 if(GetCurrentEvent()>0 && (GetCurrentEvent())%fEventsPerChannel==0)
282 fCurrentPulse = 0;
283
284
285}
286//_____________________________________________________________________
287//
288//EOF
289//