]>
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$ |
48a56166 | 20 | * Revision 1.10 2007/12/06 13:58:11 hristov |
21 | * Additional pritection. Do not delete the mapping, it is owned by another class | |
22 | * | |
e36e3bcf | 23 | * Revision 1.9 2007/12/06 02:19:51 jklay |
24 | * incorporated fitting procedure from testbeam analysis into AliRoot | |
25 | * | |
8cb998bd | 26 | * Revision 1.8 2007/12/05 02:30:51 jklay |
27 | * modification to read Altro mappings into AliEMCALRecParam and pass to AliEMCALRawUtils from AliEMCALReconstructor; add option to AliEMCALRawUtils to set old RCU format (for testbeam) or not | |
28 | * | |
feedcab9 | 29 | * Revision 1.7 2007/11/14 15:51:46 gustavo |
30 | * Take out few unnecessary prints | |
31 | * | |
6d0b6861 | 32 | * Revision 1.6 2007/11/01 01:23:51 mvl |
33 | * Removed call to SetOldRCUFormat, which is only needed for testbeam data | |
34 | * | |
f31c5945 | 35 | * Revision 1.5 2007/11/01 01:20:33 mvl |
36 | * Further improvement of peak finding; more robust fit | |
37 | * | |
e9dbb64a | 38 | * Revision 1.4 2007/10/31 17:15:24 mvl |
39 | * Fixed bug in raw data unpacking; Added pedestal to signal fit; Added logic to deal with high/low gain | |
40 | * | |
82cbdfca | 41 | * Revision 1.3 2007/09/27 08:36:46 mvl |
42 | * More robust setting of fit range in FitRawSignal (P. Hristov) | |
43 | * | |
5a056daa | 44 | * Revision 1.2 2007/09/03 20:55:35 jklay |
45 | * EMCAL e-by-e reconstruction methods from Cvetan | |
46 | * | |
c47157cd | 47 | * Revision 1.1 2007/03/17 19:56:38 mvl |
48 | * Moved signal shape routines from AliEMCAL to separate class AliEMCALRawUtils to streamline raw data reconstruction code. | |
49 | * */ | |
ee299369 | 50 | |
51 | //*-- Author: Marco van Leeuwen (LBL) | |
52 | #include "AliEMCALRawUtils.h" | |
21cad85c | 53 | |
ee299369 | 54 | #include "TF1.h" |
55 | #include "TGraph.h" | |
56 | #include "TSystem.h" | |
21cad85c | 57 | |
ee299369 | 58 | #include "AliLog.h" |
72c58de0 | 59 | #include "AliRun.h" |
ee299369 | 60 | #include "AliRunLoader.h" |
61 | #include "AliCaloAltroMapping.h" | |
62 | #include "AliAltroBuffer.h" | |
63 | #include "AliRawReader.h" | |
64 | #include "AliCaloRawStream.h" | |
65 | #include "AliDAQ.h" | |
21cad85c | 66 | |
feedcab9 | 67 | #include "AliEMCALRecParam.h" |
ee299369 | 68 | #include "AliEMCALLoader.h" |
69 | #include "AliEMCALGeometry.h" | |
70 | #include "AliEMCALDigitizer.h" | |
71 | #include "AliEMCALDigit.h" | |
20b636fc | 72 | #include "AliEMCAL.h" |
21cad85c | 73 | |
ee299369 | 74 | ClassImp(AliEMCALRawUtils) |
21cad85c | 75 | |
ee299369 | 76 | // Signal shape parameters |
21cad85c | 77 | Double_t AliEMCALRawUtils::fgTimeBinWidth = 100E-9 ; // each sample is 100 ns |
ee299369 | 78 | Double_t AliEMCALRawUtils::fgTimeTrigger = 1.5E-6 ; // 15 time bins ~ 1.5 musec |
79 | ||
80 | // some digitization constants | |
81 | Int_t AliEMCALRawUtils::fgThreshold = 1; | |
82 | Int_t AliEMCALRawUtils::fgDDLPerSuperModule = 2; // 2 ddls per SuperModule | |
fe93d365 | 83 | Int_t AliEMCALRawUtils::fgPedestalValue = 32; // pedestal value for digits2raw |
84 | Double_t AliEMCALRawUtils::fgFEENoise = 3.; // 3 ADC channels of noise (sampled) | |
ee299369 | 85 | |
8cb998bd | 86 | AliEMCALRawUtils::AliEMCALRawUtils() |
b4133f05 | 87 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), |
88 | fNPedSamples(0), fGeom(0), fOption("") | |
8cb998bd | 89 | { |
b4133f05 | 90 | |
91 | //These are default parameters. | |
92 | //Can be re-set from without with setter functions | |
ee299369 | 93 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) |
b4133f05 | 94 | fOrder = 2; // order of gamma fn |
95 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
96 | fNoiseThreshold = 3; | |
97 | fNPedSamples = 5; | |
65bdc82f | 98 | |
99 | //Get Mapping RCU files from the AliEMCALRecParam | |
100 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
101 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
102 | ||
21cad85c | 103 | for(Int_t i = 0; i < 4; i++) { |
65bdc82f | 104 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
105 | } | |
106 | ||
72c58de0 | 107 | //To make sure we match with the geometry in a simulation file, |
108 | //let's try to get it first. If not, take the default geometry | |
33c3c91a | 109 | AliRunLoader *rl = AliRunLoader::Instance(); |
c61fe3b4 | 110 | if(!rl) AliError("Cannot find RunLoader!"); |
72c58de0 | 111 | if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) { |
112 | fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); | |
113 | } else { | |
114 | AliInfo(Form("Using default geometry in raw reco")); | |
937d0661 | 115 | fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); |
65bdc82f | 116 | } |
117 | ||
72c58de0 | 118 | if(!fGeom) AliFatal(Form("Could not get geometry!")); |
119 | ||
65bdc82f | 120 | } |
121 | ||
122 | //____________________________________________________________________________ | |
5544799a | 123 | AliEMCALRawUtils::AliEMCALRawUtils(AliEMCALGeometry *pGeometry) |
124 | : fHighLowGainFactor(0.), fOrder(0), fTau(0.), fNoiseThreshold(0), | |
125 | fNPedSamples(0), fGeom(pGeometry), fOption("") | |
126 | { | |
127 | // | |
128 | // Initialize with the given geometry - constructor required by HLT | |
129 | // HLT does not use/support AliRunLoader(s) instances | |
130 | // This is a minimum intervention solution | |
131 | // Comment by MPloskon@lbl.gov | |
132 | // | |
133 | ||
134 | //These are default parameters. | |
135 | //Can be re-set from without with setter functions | |
136 | fHighLowGainFactor = 16. ; // adjusted for a low gain range of 82 GeV (10 bits) | |
137 | fOrder = 2; // order of gamma fn | |
138 | fTau = 2.35; // in units of timebin, from CERN 2007 testbeam | |
139 | fNoiseThreshold = 3; | |
140 | fNPedSamples = 5; | |
141 | ||
142 | //Get Mapping RCU files from the AliEMCALRecParam | |
143 | const TObjArray* maps = AliEMCALRecParam::GetMappings(); | |
144 | if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!"); | |
145 | ||
21cad85c | 146 | for(Int_t i = 0; i < 4; i++) { |
5544799a | 147 | fMapping[i] = (AliAltroMapping*)maps->At(i); |
148 | } | |
149 | ||
150 | if(!fGeom) AliFatal(Form("Could not get geometry!")); | |
151 | ||
152 | } | |
153 | ||
154 | //____________________________________________________________________________ | |
65bdc82f | 155 | AliEMCALRawUtils::AliEMCALRawUtils(const AliEMCALRawUtils& rawU) |
156 | : TObject(), | |
157 | fHighLowGainFactor(rawU.fHighLowGainFactor), | |
b4133f05 | 158 | fOrder(rawU.fOrder), |
159 | fTau(rawU.fTau), | |
160 | fNoiseThreshold(rawU.fNoiseThreshold), | |
161 | fNPedSamples(rawU.fNPedSamples), | |
65bdc82f | 162 | fGeom(rawU.fGeom), |
163 | fOption(rawU.fOption) | |
164 | { | |
165 | //copy ctor | |
166 | fMapping[0] = rawU.fMapping[0]; | |
167 | fMapping[1] = rawU.fMapping[1]; | |
21cad85c | 168 | fMapping[2] = rawU.fMapping[2]; |
169 | fMapping[3] = rawU.fMapping[3]; | |
65bdc82f | 170 | } |
171 | ||
172 | //____________________________________________________________________________ | |
173 | AliEMCALRawUtils& AliEMCALRawUtils::operator =(const AliEMCALRawUtils &rawU) | |
174 | { | |
175 | //assignment operator | |
176 | ||
177 | if(this != &rawU) { | |
178 | fHighLowGainFactor = rawU.fHighLowGainFactor; | |
b4133f05 | 179 | fOrder = rawU.fOrder; |
180 | fTau = rawU.fTau; | |
181 | fNoiseThreshold = rawU.fNoiseThreshold; | |
182 | fNPedSamples = rawU.fNPedSamples; | |
65bdc82f | 183 | fGeom = rawU.fGeom; |
184 | fOption = rawU.fOption; | |
185 | fMapping[0] = rawU.fMapping[0]; | |
186 | fMapping[1] = rawU.fMapping[1]; | |
21cad85c | 187 | fMapping[2] = rawU.fMapping[2]; |
188 | fMapping[3] = rawU.fMapping[3]; | |
65bdc82f | 189 | } |
190 | ||
191 | return *this; | |
192 | ||
ee299369 | 193 | } |
65bdc82f | 194 | |
ee299369 | 195 | //____________________________________________________________________________ |
196 | AliEMCALRawUtils::~AliEMCALRawUtils() { | |
65bdc82f | 197 | |
ee299369 | 198 | } |
65bdc82f | 199 | |
ee299369 | 200 | //____________________________________________________________________________ |
65bdc82f | 201 | void AliEMCALRawUtils::Digits2Raw() |
ee299369 | 202 | { |
203 | // convert digits of the current event to raw data | |
204 | ||
33c3c91a | 205 | AliRunLoader *rl = AliRunLoader::Instance(); |
ee299369 | 206 | AliEMCALLoader *loader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL")); |
207 | ||
208 | // get the digits | |
209 | loader->LoadDigits("EMCAL"); | |
210 | loader->GetEvent(); | |
211 | TClonesArray* digits = loader->Digits() ; | |
212 | ||
213 | if (!digits) { | |
214 | Warning("Digits2Raw", "no digits found !"); | |
215 | return; | |
216 | } | |
65bdc82f | 217 | |
ee299369 | 218 | static const Int_t nDDL = 12*2; // 12 SM hardcoded for now. Buffers allocated dynamically, when needed, so just need an upper limit here |
219 | AliAltroBuffer* buffers[nDDL]; | |
220 | for (Int_t i=0; i < nDDL; i++) | |
221 | buffers[i] = 0; | |
222 | ||
223 | Int_t adcValuesLow[fgkTimeBins]; | |
224 | Int_t adcValuesHigh[fgkTimeBins]; | |
225 | ||
ee299369 | 226 | // loop over digits (assume ordered digits) |
227 | for (Int_t iDigit = 0; iDigit < digits->GetEntries(); iDigit++) { | |
228 | AliEMCALDigit* digit = dynamic_cast<AliEMCALDigit *>(digits->At(iDigit)) ; | |
229 | if (digit->GetAmp() < fgThreshold) | |
230 | continue; | |
231 | ||
232 | //get cell indices | |
233 | Int_t nSM = 0; | |
234 | Int_t nIphi = 0; | |
235 | Int_t nIeta = 0; | |
236 | Int_t iphi = 0; | |
237 | Int_t ieta = 0; | |
238 | Int_t nModule = 0; | |
65bdc82f | 239 | fGeom->GetCellIndex(digit->GetId(), nSM, nModule, nIphi, nIeta); |
240 | fGeom->GetCellPhiEtaIndexInSModule(nSM, nModule, nIphi, nIeta,iphi, ieta) ; | |
ee299369 | 241 | |
21cad85c | 242 | //Check which is the RCU, 0 or 1, of the cell. |
ee299369 | 243 | Int_t iRCU = -111; |
244 | //RCU0 | |
245 | if (0<=iphi&&iphi<8) iRCU=0; // first cable row | |
246 | else if (8<=iphi&&iphi<16 && 0<=ieta&&ieta<24) iRCU=0; // first half; | |
247 | //second cable row | |
248 | //RCU1 | |
249 | else if(8<=iphi&&iphi<16 && 24<=ieta&&ieta<48) iRCU=1; // second half; | |
250 | //second cable row | |
251 | else if(16<=iphi&&iphi<24) iRCU=1; // third cable row | |
21cad85c | 252 | |
253 | if (nSM%2==1) iRCU = 1 - iRCU; // swap for odd=C side, to allow us to cable both sides the same | |
254 | ||
e36e3bcf | 255 | if (iRCU<0) |
256 | Fatal("Digits2Raw()","Non-existent RCU number: %d", iRCU); | |
ee299369 | 257 | |
258 | //Which DDL? | |
259 | Int_t iDDL = fgDDLPerSuperModule* nSM + iRCU; | |
260 | if (iDDL >= nDDL) | |
261 | Fatal("Digits2Raw()","Non-existent DDL board number: %d", iDDL); | |
262 | ||
263 | if (buffers[iDDL] == 0) { | |
264 | // open new file and write dummy header | |
265 | TString fileName = AliDAQ::DdlFileName("EMCAL",iDDL); | |
21cad85c | 266 | //Select mapping file RCU0A, RCU0C, RCU1A, RCU1C |
267 | Int_t iRCUside=iRCU+(nSM%2)*2; | |
268 | //iRCU=0 and even (0) SM -> RCU0A.data 0 | |
269 | //iRCU=1 and even (0) SM -> RCU1A.data 1 | |
270 | //iRCU=0 and odd (1) SM -> RCU0C.data 2 | |
271 | //iRCU=1 and odd (1) SM -> RCU1C.data 3 | |
272 | //cout<<" nSM "<<nSM<<"; iRCU "<<iRCU<<"; iRCUside "<<iRCUside<<endl; | |
273 | buffers[iDDL] = new AliAltroBuffer(fileName.Data(),fMapping[iRCUside]); | |
ee299369 | 274 | buffers[iDDL]->WriteDataHeader(kTRUE, kFALSE); //Dummy; |
275 | } | |
276 | ||
277 | // out of time range signal (?) | |
278 | if (digit->GetTimeR() > GetRawFormatTimeMax() ) { | |
279 | AliInfo("Signal is out of time range.\n"); | |
280 | buffers[iDDL]->FillBuffer((Int_t)digit->GetAmp()); | |
281 | buffers[iDDL]->FillBuffer(GetRawFormatTimeBins() ); // time bin | |
282 | buffers[iDDL]->FillBuffer(3); // bunch length | |
283 | buffers[iDDL]->WriteTrailer(3, ieta, iphi, nSM); // trailer | |
284 | // calculate the time response function | |
285 | } else { | |
286 | Bool_t lowgain = RawSampledResponse(digit->GetTimeR(), digit->GetAmp(), adcValuesHigh, adcValuesLow) ; | |
287 | if (lowgain) | |
288 | buffers[iDDL]->WriteChannel(ieta, iphi, 0, GetRawFormatTimeBins(), adcValuesLow, fgThreshold); | |
289 | else | |
290 | buffers[iDDL]->WriteChannel(ieta,iphi, 1, GetRawFormatTimeBins(), adcValuesHigh, fgThreshold); | |
291 | } | |
292 | } | |
293 | ||
294 | // write headers and close files | |
295 | for (Int_t i=0; i < nDDL; i++) { | |
296 | if (buffers[i]) { | |
297 | buffers[i]->Flush(); | |
298 | buffers[i]->WriteDataHeader(kFALSE, kFALSE); | |
299 | delete buffers[i]; | |
300 | } | |
301 | } | |
65bdc82f | 302 | |
ee299369 | 303 | loader->UnloadDigits(); |
304 | } | |
305 | ||
306 | //____________________________________________________________________________ | |
65bdc82f | 307 | void AliEMCALRawUtils::Raw2Digits(AliRawReader* reader,TClonesArray *digitsArr) |
ee299369 | 308 | { |
65bdc82f | 309 | // convert raw data of the current event to digits |
ee299369 | 310 | |
c47157cd | 311 | digitsArr->Clear(); |
ee299369 | 312 | |
c47157cd | 313 | if (!digitsArr) { |
ee299369 | 314 | Error("Raw2Digits", "no digits found !"); |
315 | return; | |
316 | } | |
317 | if (!reader) { | |
318 | Error("Raw2Digits", "no raw reader found !"); | |
319 | return; | |
320 | } | |
321 | ||
65bdc82f | 322 | AliCaloRawStream in(reader,"EMCAL",fMapping); |
ee299369 | 323 | // Select EMCAL DDL's; |
324 | reader->Select("EMCAL"); | |
feedcab9 | 325 | |
8cb998bd | 326 | //Updated fitting routine from 2007 beam test takes into account |
327 | //possibility of two peaks in data and selects first one for fitting | |
328 | //Also sets some of the starting parameters based on the shape of the | |
329 | //given raw signal being fit | |
330 | ||
331 | TF1 * signalF = new TF1("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); | |
b4133f05 | 332 | signalF->SetParameters(10.,0.,fTau,fOrder,5.); //set all defaults once, just to be safe |
8cb998bd | 333 | signalF->SetParNames("amp","t0","tau","N","ped"); |
b4133f05 | 334 | signalF->SetParameter(2,fTau); // tau in units of time bin |
8cb998bd | 335 | signalF->SetParLimits(2,2,-1); |
b4133f05 | 336 | signalF->SetParameter(3,fOrder); // order |
8cb998bd | 337 | signalF->SetParLimits(3,2,-1); |
48a56166 | 338 | |
ee299369 | 339 | Int_t id = -1; |
82cbdfca | 340 | Float_t time = 0. ; |
341 | Float_t amp = 0. ; | |
ee299369 | 342 | |
8cb998bd | 343 | //Graph to hold data we will fit (should be converted to an array |
344 | //later to speed up processing | |
345 | TGraph * gSig = new TGraph(GetRawFormatTimeBins()); | |
ee299369 | 346 | |
82cbdfca | 347 | Int_t readOk = 1; |
ee299369 | 348 | Int_t lowGain = 0; |
349 | ||
82cbdfca | 350 | while (readOk && in.GetModule() < 0) |
351 | readOk = in.Next(); // Go to first digit | |
e9dbb64a | 352 | |
8cb998bd | 353 | Int_t col = 0; |
354 | Int_t row = 0; | |
355 | ||
82cbdfca | 356 | while (readOk) { |
21cad85c | 357 | |
65bdc82f | 358 | id = fGeom->GetAbsCellIdFromCellIndexes(in.GetModule(), in.GetRow(), in.GetColumn()) ; |
ee299369 | 359 | lowGain = in.IsLowGain(); |
82cbdfca | 360 | Int_t maxTime = in.GetTime(); // timebins come in reverse order |
e9dbb64a | 361 | if (maxTime < 0 || maxTime >= GetRawFormatTimeBins()) { |
362 | AliWarning(Form("Invalid time bin %d",maxTime)); | |
363 | maxTime = GetRawFormatTimeBins(); | |
364 | } | |
82cbdfca | 365 | gSig->Set(maxTime+1); |
366 | // There is some kind of zero-suppression in the raw data, | |
367 | // so set up the TGraph in advance | |
368 | for (Int_t i=0; i < maxTime; i++) { | |
8cb998bd | 369 | gSig->SetPoint(i, i , 0); |
82cbdfca | 370 | } |
ee299369 | 371 | |
82cbdfca | 372 | Int_t iTime = 0; |
ee299369 | 373 | do { |
82cbdfca | 374 | if (in.GetTime() >= gSig->GetN()) { |
375 | AliWarning("Too many time bins"); | |
376 | gSig->Set(in.GetTime()); | |
ee299369 | 377 | } |
8cb998bd | 378 | col = in.GetColumn(); |
379 | row = in.GetRow(); | |
380 | ||
381 | gSig->SetPoint(in.GetTime(), in.GetTime(), in.GetSignal()) ; | |
fb070798 | 382 | |
82cbdfca | 383 | if (in.GetTime() > maxTime) |
384 | maxTime = in.GetTime(); | |
ee299369 | 385 | iTime++; |
82cbdfca | 386 | } while ((readOk = in.Next()) && !in.IsNewHWAddress()); |
ee299369 | 387 | |
388 | FitRaw(gSig, signalF, amp, time) ; | |
ee299369 | 389 | |
c61fe3b4 | 390 | if (amp > 0 && amp < 2000) { //check both high and low end of |
391 | //result, 2000 is somewhat arbitrary | |
82cbdfca | 392 | AliDebug(2,Form("id %d lowGain %d amp %g", id, lowGain, amp)); |
fb070798 | 393 | |
82cbdfca | 394 | AddDigit(digitsArr, id, lowGain, (Int_t)amp, time); |
ee299369 | 395 | } |
ee299369 | 396 | |
397 | // Reset graph | |
82cbdfca | 398 | for (Int_t index = 0; index < gSig->GetN(); index++) { |
8cb998bd | 399 | gSig->SetPoint(index, index, 0) ; |
ee299369 | 400 | } |
48a56166 | 401 | // Reset starting parameters for fit function |
b4133f05 | 402 | signalF->SetParameters(10.,0.,fTau,fOrder,5.); //reset all defaults just to be safe |
48a56166 | 403 | |
82cbdfca | 404 | }; // EMCAL entries loop |
ee299369 | 405 | |
406 | delete signalF ; | |
407 | delete gSig; | |
408 | ||
409 | return ; | |
410 | } | |
411 | ||
412 | //____________________________________________________________________________ | |
82cbdfca | 413 | void AliEMCALRawUtils::AddDigit(TClonesArray *digitsArr, Int_t id, Int_t lowGain, Int_t amp, Float_t time) { |
414 | // | |
415 | // Add a new digit. | |
416 | // This routine checks whether a digit exists already for this tower | |
417 | // and then decides whether to use the high or low gain info | |
418 | // | |
419 | // Called by Raw2Digits | |
420 | ||
421 | AliEMCALDigit *digit = 0, *tmpdigit = 0; | |
422 | ||
423 | TIter nextdigit(digitsArr); | |
424 | while (digit == 0 && (tmpdigit = (AliEMCALDigit*) nextdigit())) { | |
425 | if (tmpdigit->GetId() == id) | |
426 | digit = tmpdigit; | |
427 | } | |
428 | ||
429 | if (!digit) { // no digit existed for this tower; create one | |
430 | if (lowGain) | |
431 | amp = Int_t(fHighLowGainFactor * amp); | |
432 | Int_t idigit = digitsArr->GetEntries(); | |
433 | new((*digitsArr)[idigit]) AliEMCALDigit( -1, -1, id, amp, time, idigit) ; | |
434 | } | |
435 | else { // a digit already exists, check range | |
b4133f05 | 436 | // (use high gain if signal < cut value, otherwise low gain) |
82cbdfca | 437 | if (lowGain) { // new digit is low gain |
b4133f05 | 438 | if (digit->GetAmp() > fgkOverflowCut) { // use if stored digit is out of range |
82cbdfca | 439 | digit->SetAmp(Int_t(fHighLowGainFactor * amp)); |
440 | digit->SetTime(time); | |
441 | } | |
442 | } | |
b4133f05 | 443 | else if (amp < fgkOverflowCut) { // new digit is high gain; use if not out of range |
82cbdfca | 444 | digit->SetAmp(amp); |
445 | digit->SetTime(time); | |
446 | } | |
447 | } | |
448 | } | |
449 | ||
450 | //____________________________________________________________________________ | |
451 | void AliEMCALRawUtils::FitRaw(TGraph * gSig, TF1* signalF, Float_t & amp, Float_t & time) | |
ee299369 | 452 | { |
453 | // Fits the raw signal time distribution; from AliEMCALGetter | |
454 | ||
ee299369 | 455 | amp = time = 0. ; |
82cbdfca | 456 | Double_t ped = 0; |
457 | Int_t nPed = 0; | |
ee299369 | 458 | |
b4133f05 | 459 | for (Int_t index = 0; index < fNPedSamples; index++) { |
e9dbb64a | 460 | Double_t ttime, signal; |
82cbdfca | 461 | gSig->GetPoint(index, ttime, signal) ; |
462 | if (signal > 0) { | |
463 | ped += signal; | |
464 | nPed++; | |
ee299369 | 465 | } |
466 | } | |
e9dbb64a | 467 | |
82cbdfca | 468 | if (nPed > 0) |
469 | ped /= nPed; | |
e9dbb64a | 470 | else { |
8cb998bd | 471 | AliWarning("Could not determine pedestal"); |
82cbdfca | 472 | ped = 10; // put some small value as first guess |
e9dbb64a | 473 | } |
82cbdfca | 474 | |
e9dbb64a | 475 | Int_t max_found = 0; |
476 | Int_t i_max = 0; | |
477 | Float_t max = -1; | |
8cb998bd | 478 | Float_t max_fit = gSig->GetN(); |
e9dbb64a | 479 | Float_t min_after_sig = 9999; |
8cb998bd | 480 | Int_t tmin_after_sig = gSig->GetN(); |
e9dbb64a | 481 | Int_t n_ped_after_sig = 0; |
fb070798 | 482 | Int_t plateau_width = 0; |
483 | Int_t plateau_start = 9999; | |
c61fe3b4 | 484 | Float_t Cut = 0.3; |
e9dbb64a | 485 | |
b4133f05 | 486 | for (Int_t i=fNPedSamples; i < gSig->GetN(); i++) { |
e9dbb64a | 487 | Double_t ttime, signal; |
488 | gSig->GetPoint(i, ttime, signal) ; | |
489 | if (!max_found && signal > max) { | |
490 | i_max = i; | |
e9dbb64a | 491 | max = signal; |
492 | } | |
b4133f05 | 493 | else if ( max > ped + fNoiseThreshold ) { |
e9dbb64a | 494 | max_found = 1; |
495 | min_after_sig = signal; | |
8cb998bd | 496 | tmin_after_sig = i; |
e9dbb64a | 497 | } |
498 | if (max_found) { | |
499 | if ( signal < min_after_sig) { | |
500 | min_after_sig = signal; | |
8cb998bd | 501 | tmin_after_sig = i; |
e9dbb64a | 502 | } |
503 | if (i > tmin_after_sig + 5) { // Two close peaks; end fit at minimum | |
504 | max_fit = tmin_after_sig; | |
505 | break; | |
506 | } | |
c61fe3b4 | 507 | if ( signal < Cut*max){ //stop fit at 30% amplitude(avoid the pulse shape falling edge) |
508 | max_fit = i; | |
509 | break; | |
510 | } | |
b4133f05 | 511 | if ( signal < ped + fNoiseThreshold) |
e9dbb64a | 512 | n_ped_after_sig++; |
513 | if (n_ped_after_sig >= 5) { // include 5 pedestal bins after peak | |
8cb998bd | 514 | max_fit = i; |
e9dbb64a | 515 | break; |
516 | } | |
5a056daa | 517 | } |
fb070798 | 518 | //Add check on plateau |
519 | if (signal >= fgkRawSignalOverflow - fNoiseThreshold) { | |
520 | if(plateau_width == 0) plateau_start = i; | |
521 | plateau_width++; | |
522 | } | |
523 | } | |
524 | ||
525 | if(plateau_width > 0) { | |
c61fe3b4 | 526 | for(int j = 0; j < plateau_width; j++) { |
fb070798 | 527 | //Note, have to remove the same point N times because after each |
528 | //remove, the positions of all subsequent points have shifted down | |
c61fe3b4 | 529 | gSig->RemovePoint(plateau_start); |
fb070798 | 530 | } |
82cbdfca | 531 | } |
5a056daa | 532 | |
b4133f05 | 533 | if ( max - ped > fNoiseThreshold ) { // else its noise |
e9dbb64a | 534 | AliDebug(2,Form("Fitting max %d ped %d", max, ped)); |
8cb998bd | 535 | signalF->SetRange(0,max_fit); |
536 | ||
537 | if(max-ped > 50) | |
538 | signalF->SetParLimits(2,1,3); | |
539 | ||
540 | signalF->SetParameter(4, ped) ; | |
541 | signalF->SetParameter(1, i_max); | |
542 | signalF->SetParameter(0, max); | |
543 | ||
544 | gSig->Fit(signalF, "QROW"); // Note option 'W': equal errors on all points | |
545 | amp = signalF->GetParameter(0); | |
546 | time = signalF->GetParameter(1)*GetRawFormatTimeBinWidth() - fgTimeTrigger; | |
ee299369 | 547 | } |
548 | return; | |
549 | } | |
550 | //__________________________________________________________________ | |
551 | Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par) | |
552 | { | |
8cb998bd | 553 | // Matches version used in 2007 beam test |
554 | // | |
ee299369 | 555 | // Shape of the electronics raw reponse: |
556 | // It is a semi-gaussian, 2nd order Gamma function of the general form | |
557 | // | |
558 | // t' = (t - t0 + tau) / tau | |
559 | // F = A * t**N * exp( N * ( 1 - t) ) for t >= 0 | |
560 | // F = 0 for t < 0 | |
561 | // | |
562 | // parameters: | |
8cb998bd | 563 | // A: par[0] // Amplitude = peak value |
564 | // t0: par[1] | |
565 | // tau: par[2] | |
566 | // N: par[3] | |
567 | // ped: par[4] | |
ee299369 | 568 | // |
569 | Double_t signal ; | |
8cb998bd | 570 | Double_t tau =par[2]; |
571 | Double_t N =par[3]; | |
572 | Double_t ped = par[4]; | |
573 | Double_t xx = ( x[0] - par[1] + tau ) / tau ; | |
ee299369 | 574 | |
5a056daa | 575 | if (xx <= 0) |
8cb998bd | 576 | signal = ped ; |
ee299369 | 577 | else { |
8cb998bd | 578 | signal = ped + par[0] * TMath::Power(xx , N) * TMath::Exp(N * (1 - xx )) ; |
ee299369 | 579 | } |
580 | return signal ; | |
581 | } | |
582 | ||
583 | //__________________________________________________________________ | |
584 | Bool_t AliEMCALRawUtils::RawSampledResponse( | |
585 | const Double_t dtime, const Double_t damp, Int_t * adcH, Int_t * adcL) const | |
586 | { | |
587 | // for a start time dtime and an amplitude damp given by digit, | |
588 | // calculates the raw sampled response AliEMCAL::RawResponseFunction | |
589 | ||
ee299369 | 590 | Bool_t lowGain = kFALSE ; |
591 | ||
48a56166 | 592 | // A: par[0] // Amplitude = peak value |
593 | // t0: par[1] | |
594 | // tau: par[2] | |
595 | // N: par[3] | |
596 | // ped: par[4] | |
597 | ||
56e13066 | 598 | TF1 signalF("signal", RawResponseFunction, 0, GetRawFormatTimeBins(), 5); |
48a56166 | 599 | signalF.SetParameter(0, damp) ; |
56e13066 | 600 | signalF.SetParameter(1, (dtime + fgTimeTrigger)/fgTimeBinWidth) ; |
b4133f05 | 601 | signalF.SetParameter(2, fTau) ; |
602 | signalF.SetParameter(3, fOrder); | |
fe93d365 | 603 | signalF.SetParameter(4, fgPedestalValue); |
ee299369 | 604 | |
605 | for (Int_t iTime = 0; iTime < GetRawFormatTimeBins(); iTime++) { | |
56e13066 | 606 | Double_t signal = signalF.Eval(iTime) ; |
fe93d365 | 607 | |
608 | //According to Terry Awes, 13-Apr-2008 | |
609 | //add gaussian noise in quadrature to each sample | |
610 | //Double_t noise = gRandom->Gaus(0.,fgFEENoise); | |
611 | //signal = sqrt(signal*signal + noise*noise); | |
612 | ||
ee299369 | 613 | adcH[iTime] = static_cast<Int_t>(signal + 0.5) ; |
b4133f05 | 614 | if ( adcH[iTime] > fgkRawSignalOverflow ){ // larger than 10 bits |
615 | adcH[iTime] = fgkRawSignalOverflow ; | |
ee299369 | 616 | lowGain = kTRUE ; |
617 | } | |
618 | ||
619 | signal /= fHighLowGainFactor; | |
620 | ||
621 | adcL[iTime] = static_cast<Int_t>(signal + 0.5) ; | |
b4133f05 | 622 | if ( adcL[iTime] > fgkRawSignalOverflow) // larger than 10 bits |
623 | adcL[iTime] = fgkRawSignalOverflow ; | |
ee299369 | 624 | } |
625 | return lowGain ; | |
626 | } |