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