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