]>
Commit | Line | Data |
---|---|---|
ee299369 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, 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 | /* $Id$ */ | |
17 | /* History of cvs commits: | |
18 | * | |
19 | * $Log$ */ | |
20 | ||
21 | //*-- Author: Marco van Leeuwen (LBL) | |
22 | #include "AliEMCALRawUtils.h" | |
23 | ||
24 | #include "TF1.h" | |
25 | #include "TGraph.h" | |
26 | #include "TSystem.h" | |
27 | ||
28 | #include "AliLog.h" | |
29 | #include "AliRunLoader.h" | |
30 | #include "AliCaloAltroMapping.h" | |
31 | #include "AliAltroBuffer.h" | |
32 | #include "AliRawReader.h" | |
33 | #include "AliCaloRawStream.h" | |
34 | #include "AliDAQ.h" | |
35 | ||
36 | #include "AliEMCALLoader.h" | |
37 | #include "AliEMCALGeometry.h" | |
38 | #include "AliEMCALDigitizer.h" | |
39 | #include "AliEMCALDigit.h" | |
40 | ||
41 | ||
42 | ClassImp(AliEMCALRawUtils) | |
43 | ||
44 | // Signal shape parameters | |
45 | Int_t AliEMCALRawUtils::fgOrder = 2 ; // Order of gamma function | |
46 | Double_t AliEMCALRawUtils::fgTimeMax = 2.56E-5 ; // each sample is over 100 ns fTimeMax/fTimeBins | |
47 | Double_t AliEMCALRawUtils::fgTau = 165E-9 ; // 165 ns (from testbeam; not very accurate) | |
48 | Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec | |
49 | ||
50 | // some digitization constants | |
51 | Int_t AliEMCALRawUtils::fgThreshold = 1; | |
52 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule | |
53 | ||
54 | AliEMCALRawUtils::AliEMCALRawUtils(): fHighLowGainFactor(0.) { | |
55 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) | |
56 | } | |
57 | //____________________________________________________________________________ | |
58 | AliEMCALRawUtils::~AliEMCALRawUtils() { | |
59 | } | |
60 | //____________________________________________________________________________ | |
61 | void AliEMCALRawUtils::Digits2Raw() | |
62 | { | |
63 | // convert digits of the current event to raw data | |
64 | ||
65 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); | |
66 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); | |
67 | ||
68 | // get the digits | |
69 | loader->LoadDigits("EMCAL"); | |
70 | loader->GetEvent(); | |
71 | TClonesArray* digits = loader->Digits() ; | |
72 | ||
73 | if (!digits) { | |
74 | Warning("Digits2Raw", "no digits found !"); | |
75 | return; | |
76 | } | |
77 | ||
78 | // get the geometry | |
79 | AliEMCALGeometry* geom = AliEMCALGeometry::GetInstance(); | |
80 | if (!geom) { | |
81 | AliError(Form("No geometry found !")); | |
82 | return; | |
83 | } | |
84 | ||
85 | static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here | |
86 | AliAltroBuffer* buffers[nDDL]; | |
87 | for (Int_t i=0; i < nDDL; i++) | |
88 | buffers[i] = 0; | |
89 | ||
90 | Int_t adcValuesLow[fgkTimeBins]; | |
91 | Int_t adcValuesHigh[fgkTimeBins]; | |
92 | ||
93 | //Load Mapping RCU files once | |
94 | TString path = gSystem->Getenv("ALICE_ROOT"); | |
95 | path += "/EMCAL/mapping/RCU"; | |
96 | TString path0 = path+"0.data";//This file will change in future | |
97 | TString path1 = path+"1.data";//This file will change in future | |
98 | AliAltroMapping * mapping[2] ; // For the moment only 2 | |
99 | mapping[0] = new AliCaloAltroMapping(path0.Data()); | |
100 | mapping[1] = new AliCaloAltroMapping(path1.Data()); | |
101 | ||
102 | // loop over digits (assume ordered digits) | |
103 | for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { | |
104 | AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ; | |
105 | if (digit->GetAmp() < fgThreshold) | |
106 | continue; | |
107 | ||
108 | //get cell indices | |
109 | Int_t nSM = 0; | |
110 | Int_t nIphi = 0; | |
111 | Int_t nIeta = 0; | |
112 | Int_t iphi = 0; | |
113 | Int_t ieta = 0; | |
114 | Int_t nModule = 0; | |
115 | geom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); | |
116 | geom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; | |
117 | ||
118 | //Check which is the RCU of the cell. | |
119 | Int_t iRCU = -111; | |
120 | //RCU0 | |
121 | if (0<=iphi&&iphi<8) iRCU=0; // first cable row | |
122 | else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; | |
123 | //second cable row | |
124 | //RCU1 | |
125 | else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; | |
126 | //second cable row | |
127 | else if(16<=iphi&&iphi<24) iRCU=1; // third cable row | |
128 | ||
129 | //Which DDL? | |
130 | Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; | |
131 | if (iDDL >= nDDL) | |
132 | Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); | |
133 | ||
134 | if (buffers[iDDL] == 0) { | |
135 | // open new file and write dummy header | |
136 | TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); | |
137 | buffers[iDDL] = new AliAltroBuffer(fileName.Data(),mapping[iRCU]); | |
138 | buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; | |
139 | } | |
140 | ||
141 | // out of time range signal (?) | |
142 | if (digit->GetTimeR() > GetRawFormatTimeMax() ) { | |
143 | AliInfo("Signal is out of time range.\n"); | |
144 | buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp()); | |
145 | buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin | |
146 | buffers[iDDL]->FillBuffer(3); // bunch length | |
147 | buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer | |
148 | // calculate the time response function | |
149 | } else { | |
150 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; | |
151 | if (lowgain) | |
152 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold); | |
153 | else | |
154 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold); | |
155 | } | |
156 | } | |
157 | ||
158 | // write headers and close files | |
159 | for (Int_t i=0; i < nDDL; i++) { | |
160 | if (buffers[i]) { | |
161 | buffers[i]->Flush(); | |
162 | buffers[i]->WriteDataHeader(kFALSE, kFALSE); | |
163 | delete buffers[i]; | |
164 | } | |
165 | } | |
166 | mapping[0]->Delete(); | |
167 | mapping[1]->Delete(); | |
168 | loader->UnloadDigits(); | |
169 | } | |
170 | ||
171 | //____________________________________________________________________________ | |
172 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader) | |
173 | { | |
174 | // convert raw data of the current event to digits | |
175 | AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); | |
176 | if (!geom) { | |
177 | AliError(Form("No geometry found !")); | |
178 | return; | |
179 | } | |
180 | ||
181 | AliRunLoader *rl = AliRunLoader::GetRunLoader(); | |
182 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); | |
183 | ||
184 | // get the digits | |
185 | TClonesArray* digits = loader->Digits() ; | |
186 | digits->Clear(); | |
187 | ||
188 | if (!digits) { | |
189 | Error("Raw2Digits", "no digits found !"); | |
190 | return; | |
191 | } | |
192 | if (!reader) { | |
193 | Error("Raw2Digits", "no raw reader found !"); | |
194 | return; | |
195 | } | |
196 | ||
197 | // Use AliAltroRawStream to read the ALTRO format. No need to | |
198 | // reinvent the wheel :-) | |
199 | AliCaloRawStream in(reader,"EMCAL"); | |
200 | // Select EMCAL DDL's; | |
201 | reader->Select("EMCAL"); | |
202 | ||
203 | // reading is from previously existing AliEMCALGetter.cxx | |
204 | // ReadRaw method | |
205 | TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4); | |
206 | signalF->SetParNames("Charge", "Gain", "Amplitude", "TimeZero"); | |
207 | ||
208 | Int_t id = -1; | |
209 | Int_t idigit = 0 ; | |
210 | Double_t time = 0. ; | |
211 | Double_t amp = 0. ; | |
212 | ||
213 | TGraph * gSig = new TGraph(GetRawFormatTimeBins()) ; | |
214 | ||
215 | Int_t eofReached = 0; | |
216 | Int_t lowGain = 0; | |
217 | ||
218 | in.Next(); // Go to first digit | |
219 | do { | |
220 | id = geom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; | |
221 | lowGain = in.IsLowGain(); | |
222 | gSig->SetPoint(in.GetTime(), | |
223 | in.GetTime()* GetRawFormatTimeMax() / GetRawFormatTimeBins(), | |
224 | in.GetSignal()) ; | |
225 | ||
226 | Int_t iTime = 1; | |
227 | do { | |
228 | if (!in.Next()) | |
229 | eofReached = 1; | |
230 | else { | |
231 | gSig->SetPoint(in.GetTime(), | |
232 | in.GetTime()* GetRawFormatTimeMax() / GetRawFormatTimeBins(), | |
233 | in.GetSignal()) ; | |
234 | } | |
235 | iTime++; | |
236 | } while (!eofReached && !in.IsNewRow() && !in.IsNewColumn() && !in.IsNewModule()); | |
237 | ||
238 | FitRaw(gSig, signalF, amp, time) ; | |
239 | if (lowGain) | |
240 | amp *= fHighLowGainFactor; | |
241 | ||
242 | if (amp > 0) { | |
243 | AliDebug(2,Form("id %d amp %g", id, amp)); | |
244 | new((*digits)[idigit]) AliEMCALDigit( -1, -1, id, (Int_t)amp, time, idigit) ; | |
245 | idigit++ ; | |
246 | } | |
247 | Int_t index ; | |
248 | ||
249 | // Reset graph | |
250 | for (index = 0; index < GetRawFormatTimeBins(); index++) { | |
251 | gSig->SetPoint(index, index * GetRawFormatTimeMax() / GetRawFormatTimeBins(), 0) ; | |
252 | } | |
253 | } while (!eofReached); // EMCAL entries loop | |
254 | ||
255 | delete signalF ; | |
256 | delete gSig; | |
257 | ||
258 | return ; | |
259 | } | |
260 | ||
261 | //____________________________________________________________________________ | |
262 | void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Double_t & amp, Double_t & time) | |
263 | { | |
264 | // Fits the raw signal time distribution; from AliEMCALGetter | |
265 | ||
266 | const Int_t kNoiseThreshold = 0 ; | |
267 | Double_t timezero1 = 0., timezero2 = 0., timemax = 0. ; | |
268 | Double_t signal = 0., signalmax = 0. ; | |
269 | amp = time = 0. ; | |
270 | ||
271 | timezero1 = signalmax = timemax = 0. ; | |
272 | timezero2 = GetRawFormatTimeMax(); | |
273 | Int_t index ; | |
274 | for (index = 0; index < GetRawFormatTimeBins(); index++) { | |
275 | gSig->GetPoint(index, time, signal) ; | |
276 | if (signal > kNoiseThreshold && timezero1 == 0.) | |
277 | timezero1 = time ; | |
278 | if (signal <= kNoiseThreshold && timezero1 > 0. && timezero2 == 0.) | |
279 | timezero2 = time ; | |
280 | if (signal > signalmax) { | |
281 | signalmax = signal ; | |
282 | timemax = time ; | |
283 | } | |
284 | } | |
285 | ||
286 | if ( timemax < GetRawFormatTimeMax() * 0.4 ) { // else its noise | |
287 | signalF->SetParameter(0, signalmax) ; | |
288 | signalF->SetParameter(1, timemax) ; | |
289 | gSig->Fit(signalF, "QRON", "", 0., timezero2); //, "QRON") ; | |
290 | amp = signalF->GetParameter(0) ; | |
291 | time = signalF->GetParameter(1) - fgTimeTrigger; | |
292 | } | |
293 | return; | |
294 | } | |
295 | //__________________________________________________________________ | |
296 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) | |
297 | { | |
298 | // Shape of the electronics raw reponse: | |
299 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
300 | // | |
301 | // t' = (t - t0 + tau) / tau | |
302 | // F = A * t**N * exp( N * ( 1 - t) ) for t >= 0 | |
303 | // F = 0 for t < 0 | |
304 | // | |
305 | // parameters: | |
306 | // A: par[0] // Amplitude = peak value | |
307 | // t0: par[1] | |
308 | // tau: fgTau | |
309 | // N: fgOrder | |
310 | // | |
311 | Double_t signal ; | |
312 | Double_t xx = ( x[0] - par[1] + fgTau ) / fgTau ; | |
313 | ||
314 | if (xx < 0) | |
315 | signal = 0. ; | |
316 | else { | |
317 | signal = par[0] * TMath::Power(xx , fgOrder) * TMath::Exp(fgOrder * (1 - xx )) ; | |
318 | } | |
319 | return signal ; | |
320 | } | |
321 | ||
322 | //__________________________________________________________________ | |
323 | Bool_t AliEMCALRawUtils::RawSampledResponse( | |
324 | const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const | |
325 | { | |
326 | // for a start time dtime and an amplitude damp given by digit, | |
327 | // calculates the raw sampled response AliEMCAL::RawResponseFunction | |
328 | ||
329 | const Int_t kRawSignalOverflow = 0x3FF ; | |
330 | Bool_t lowGain = kFALSE ; | |
331 | ||
332 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeMax(), 4); | |
333 | signalF.SetParameter(0, damp) ; | |
334 | signalF.SetParameter(1, dtime + fgTimeTrigger) ; | |
335 | ||
336 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { | |
337 | Double_t time = iTime * GetRawFormatTimeMax() / GetRawFormatTimeBins() ; | |
338 | Double_t signal = signalF.Eval(time) ; | |
339 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
340 | if ( adcH[iTime] > kRawSignalOverflow ){ // larger than 10 bits | |
341 | adcH[iTime] = kRawSignalOverflow ; | |
342 | lowGain = kTRUE ; | |
343 | } | |
344 | ||
345 | signal /= fHighLowGainFactor; | |
346 | ||
347 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
348 | if ( adcL[iTime] > kRawSignalOverflow) // larger than 10 bits | |
349 | adcL[iTime] = kRawSignalOverflow ; | |
350 | } | |
351 | return lowGain ; | |
352 | } |