]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - EMCAL/AliEMCALDigit.cxx
Script to draw occupancies
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigit.cxx
... / ...
CommitLineData
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
18//_________________________________________________________________________
19// EMCAL digit:
20// A Digit is the sum of the energy lost in an EMCAL Tower
21// It also stores information on Primary, and enterring particle
22// tracknumbers Digits are created using AliEMCALSDigitizer, followed
23// by AliEMCALDigitizer
24//
25//*-- Author: Sahal Yacoob (LBL)
26// based on : AliPHOSDigit
27//__________________________________________________________________________
28
29// --- ROOT system ---
30#include <Riostream.h>
31#include <TMath.h>
32
33// --- Standard library ---
34
35// --- AliRoot header files ---
36
37#include "AliEMCALDigit.h"
38#include "AliEMCALGeometry.h"
39
40ClassImp(AliEMCALDigit)
41
42//____________________________________________________________________________
43AliEMCALDigit::AliEMCALDigit() :
44AliDigitNew(),
45 fAmpFloat(0.),
46 fNSamples(0),
47 fSamples(0x0),
48 fNSamplesHG(0),
49 fSamplesHG(0x0),
50 fNprimary(0),
51 fNMaxPrimary(5),
52 fPrimary(0x0),
53 fDEPrimary(0x0),
54 fNiparent(0),
55 fNMaxiparent(5),
56 fIparent(0x0),
57 fDEParent(0x0),
58 fMaxIter(0),
59 fTime(0.),
60 fTimeR(0.),
61 fChi2(0.),
62 fNDF(0),
63 fDigitType(kUnknown),
64 fAmpCalib(-1)
65{
66 // default ctor
67
68 // Need to initialise for reading old files
69 fPrimary = new Int_t[fNMaxPrimary] ;
70 fDEPrimary = new Float_t[fNMaxPrimary] ;
71 fIparent = new Int_t[fNMaxiparent] ;
72 fDEParent = new Float_t[fNMaxiparent] ;
73 for ( Int_t i = 0; i < fNMaxPrimary ; i++) {
74 fPrimary[i] = -1 ;
75 fDEPrimary[i] = 0 ;
76 }
77
78 for ( Int_t i = 0; i < fNMaxiparent ; i++) {
79 fIparent[i] = -1 ;
80 fDEParent[i] = 0 ;
81 }
82}
83
84//____________________________________________________________________________
85AliEMCALDigit::AliEMCALDigit(Int_t primary, Int_t iparent, Int_t id, Float_t digEnergy, Float_t time, Int_t type, Int_t index, Float_t chi2, Int_t ndf, Float_t dE)
86 : AliDigitNew(),
87 fAmpFloat(digEnergy),
88 fNSamples(0),
89 fSamples(0x0),
90 fNSamplesHG(0),
91 fSamplesHG(0x0),
92 fNprimary(0),
93 fNMaxPrimary(25),
94 fPrimary(0x0),
95 fDEPrimary(0x0),
96 fNiparent(0),
97 fNMaxiparent(150),
98 fIparent(0x0),
99 fDEParent(0x0),
100 fMaxIter(5),
101 fTime(time),
102 fTimeR(time),
103 fChi2(chi2),
104 fNDF(ndf),
105 fDigitType(type),
106 fAmpCalib(-1)
107{
108 // ctor with all data
109
110 // data memebrs of the base class (AliNewDigit)
111 fAmp = 0;
112 fId = id ;
113 fIndexInList = index ;
114
115 // data member
116 fPrimary = new Int_t[fNMaxPrimary] ;
117 fDEPrimary = new Float_t[fNMaxPrimary] ;
118 fIparent = new Int_t[fNMaxiparent] ;
119 fDEParent = new Float_t[fNMaxiparent] ;
120 if( primary != -1){
121 fNprimary = 1 ;
122 fPrimary[0] = primary ;
123 fDEPrimary[0] = dE ;
124 fNiparent = 1 ;
125 fIparent[0] = iparent ;
126 fDEParent[0] = dE ;
127 }
128 else{ //If the contribution of this primary smaller than fDigitThreshold (AliEMCALv1)
129 fNprimary = 0 ;
130 fPrimary[0] = -1 ;
131 fDEPrimary[0] = 0 ;
132 fNiparent = 0 ;
133 fIparent[0] = -1 ;
134 fDEParent[0] = 0 ;
135 }
136 Int_t i ;
137 for ( i = 1; i < fNMaxPrimary ; i++) {
138 fPrimary[i] = -1 ;
139 fDEPrimary[i] = 0 ;
140 }
141
142 for ( i = 1; i< fNMaxiparent ; i++) {
143 fIparent[i] = -1 ;
144 fDEParent[i] = 0 ;
145 }
146}
147
148//____________________________________________________________________________
149AliEMCALDigit::AliEMCALDigit(const AliEMCALDigit & digit)
150 : AliDigitNew(digit),
151 fAmpFloat(digit.fAmpFloat),
152 fNSamples(digit.fNSamples),
153 fSamples(),
154 fNSamplesHG(digit.fNSamplesHG),
155 fSamplesHG(),
156 fNprimary(digit.fNprimary),
157 fNMaxPrimary(digit.fNMaxPrimary),
158 fPrimary(),
159 fDEPrimary(),
160 fNiparent(digit.fNiparent),
161 fNMaxiparent(digit.fNMaxiparent),
162 fIparent(),
163 fDEParent(),
164 fMaxIter(digit.fMaxIter),
165 fTime(digit.fTime),
166 fTimeR(digit.fTimeR),
167 fChi2(digit.fChi2),
168 fNDF(digit.fNDF),
169 fDigitType(digit.fDigitType),
170 fAmpCalib(digit.fAmpCalib)
171{
172 // copy ctor
173 // data memebrs of the base class (AliNewDigit)
174 fAmp = digit.fAmp ;
175 fId = digit.fId;
176 fIndexInList = digit.fIndexInList ;
177
178 // data members
179 if (fSamples ) delete [] fSamples ; fSamples = NULL ;
180 if (fSamplesHG ) delete [] fSamplesHG ; fSamplesHG = NULL ;
181 if (fPrimary ) delete [] fPrimary ; fPrimary = NULL ;
182 if (fDEPrimary) delete [] fDEPrimary ; fDEPrimary = NULL ;
183 if (fIparent ) delete [] fIparent ; fIparent = NULL ;
184 if (fDEParent) delete [] fDEParent ; fDEParent = NULL ;
185
186 if (fNSamples){
187 fSamples = new Int_t[fNSamples];
188 for (Int_t i=0; i < digit.fNSamples; i++) fSamples[i] = digit.fSamples[i];
189 }
190
191 if (fNSamplesHG){
192 fSamplesHG = new Int_t[fNSamplesHG];
193 for (Int_t i=0; i < digit.fNSamplesHG; i++) fSamplesHG[i] = digit.fSamplesHG[i];
194 }
195
196
197 if (fNMaxPrimary){
198 fPrimary = new Int_t [fNMaxPrimary] ;
199 fDEPrimary = new Float_t[fNMaxPrimary] ;
200 for ( Int_t i = 0; i < fNMaxPrimary ; i++) {
201 fPrimary[i] = digit.fPrimary[i] ;
202 fDEPrimary[i] = digit.fDEPrimary[i] ;
203 }
204 }
205
206 if (fNMaxiparent){
207 fIparent = new Int_t [fNMaxiparent] ;
208 fDEParent = new Float_t[fNMaxiparent] ;
209 for (Int_t j = 0; j< fNMaxiparent ; j++) {
210 fIparent[j] = digit.fIparent[j] ;
211 fDEParent[j] = digit.fDEParent[j] ;
212 }
213 }
214
215}
216
217//____________________________________________________________________________
218AliEMCALDigit::~AliEMCALDigit()
219{
220 // Delete array of primaries if any
221
222 if (fSamples ) delete [] fSamples ; fSamples = NULL ;
223 if (fSamplesHG ) delete [] fSamplesHG ; fSamplesHG = NULL ;
224 if (fPrimary ) delete [] fPrimary ; fPrimary = NULL ;
225 if (fDEPrimary) delete [] fDEPrimary ; fDEPrimary = NULL ;
226 if (fIparent ) delete [] fIparent ; fIparent = NULL ;
227 if (fDEParent) delete [] fDEParent ; fDEParent = NULL ;
228
229}
230
231
232//____________________________________________________________________________
233void AliEMCALDigit::Clear(const Option_t*)
234{
235 // Delete array of primaries if any
236 if (fSamples ) delete [] fSamples ; fSamples = NULL ;
237 if (fSamplesHG ) delete [] fSamplesHG ; fSamplesHG = NULL ;
238 if (fPrimary ) delete [] fPrimary ; fPrimary = NULL ;
239 if (fDEPrimary) delete [] fDEPrimary ; fDEPrimary = NULL ;
240 if (fIparent ) delete [] fIparent ; fIparent = NULL ;
241 if (fDEParent) delete [] fDEParent ; fDEParent = NULL ;
242
243}
244
245//____________________________________________________________________________
246Int_t AliEMCALDigit::Compare(const TObject * obj) const
247{
248 // Compares two digits with respect to its Id
249 // to sort according increasing Id
250
251 Int_t rv = 2 ;
252
253 AliEMCALDigit * digit = (AliEMCALDigit *)obj ;
254
255 Int_t iddiff = fId - digit->GetId() ;
256
257 if ( iddiff > 0 )
258 rv = 1 ;
259 else if ( iddiff < 0 )
260 rv = -1 ;
261 else
262 rv = 0 ;
263
264 return rv ;
265
266}
267
268//____________________________________________________________________________
269Float_t AliEMCALDigit::GetEta() const
270{
271 //return pseudorapidity for this digit
272 // should be change in EMCALGeometry - 19-nov-04
273 Float_t eta=-10., phi=-10.;
274 Int_t id = GetId();
275 const AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
276 g->EtaPhiFromIndex(id,eta,phi);
277 return eta ;
278}
279
280//____________________________________________________________________________
281Float_t AliEMCALDigit::GetPhi() const
282{
283 //return phi coordinate of digit
284 // should be change in EMCALGeometry - 19-nov-04
285 Float_t eta=-10., phi=-10.;
286 Int_t id = GetId();
287 const AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
288 g->EtaPhiFromIndex(id,eta,phi);
289 return phi ;
290}
291
292//____________________________________________________________________________
293Bool_t AliEMCALDigit::GetFALTROSample(const Int_t iSample, Int_t& timeBin, Int_t& amp) const
294{
295 //Get FALTRO sample in time bin iSample
296 if (iSample >= fNSamples || iSample < 0 || fDigitType==kTrigger) return kFALSE;
297
298 amp = fSamples[iSample] & 0xFFF;
299 timeBin = (fSamples[iSample] >> 12) & 0xFF;
300
301 return kTRUE;
302}
303//____________________________________________________________________________
304Bool_t AliEMCALDigit::GetALTROSampleLG(const Int_t iSample, Int_t& timeBin, Int_t& amp) const
305{
306 //Get Low Gain ALTRO sample in time bin iSample
307 if (iSample >= fNSamples || iSample < 0 || fDigitType==kLG) return kFALSE;
308
309 amp = fSamples[iSample] & 0xFFF;
310 timeBin = (fSamples[iSample] >> 12) & 0xFF;
311
312 return kTRUE;
313}
314
315//____________________________________________________________________________
316void AliEMCALDigit::SetALTROSamplesLG(const Int_t nSamples, Int_t *samples)
317{
318 //Set array of ALTRO samples, Low Gain or FALTRO
319 fNSamples = nSamples;
320 fSamples = new Int_t[fNSamples];
321 for (Int_t i=0; i < fNSamples; i++) fSamples[i] = samples[i];
322}
323
324//____________________________________________________________________________
325void AliEMCALDigit::SetALTROSamplesHG(const Int_t nSamples, Int_t *samples)
326{
327 //Set array of ALTRO samples, High Gain.
328 fNSamplesHG = nSamples;
329 fSamplesHG = new Int_t[fNSamplesHG];
330 for (Int_t i=0; i < fNSamplesHG; i++) fSamplesHG[i] = samples[i];
331}
332
333
334//____________________________________________________________________________
335Bool_t AliEMCALDigit::GetALTROSampleHG(const Int_t iSample, Int_t& timeBin, Int_t& amp) const
336{
337 //Get High Gain ALTRO sample in time bin iSample
338 if (iSample >= fNSamplesHG || iSample < 0 || fDigitType==kHG) return kFALSE;
339
340 amp = fSamplesHG[iSample] & 0xFFF;
341 timeBin = (fSamplesHG[iSample] >> 12) & 0xFF;
342
343 return kTRUE;
344}
345
346
347//____________________________________________________________________________
348Int_t AliEMCALDigit::GetPrimary(Int_t index) const
349{
350 // retrieves the primary particle number given its index in the list
351 if ( (index <= fNprimary) && (index > 0)){
352 return fPrimary[index-1] ;
353 }
354
355 return -1 ;
356}
357
358//____________________________________________________________________________
359Float_t AliEMCALDigit::GetDEPrimary(Int_t index) const
360{
361 // retrieves the primary particle energy contribution
362 // given its index in the list
363 if ( (index <= fNprimary) && (index > 0)){
364 return fDEPrimary[index-1] ;
365 }
366
367 return 0 ;
368
369}
370
371//____________________________________________________________________________
372Int_t AliEMCALDigit::GetIparent(Int_t index) const
373{
374 // retrieves the primary particle number given its index in the list
375 if ( index <= fNiparent && index > 0){
376 return fIparent[index-1] ;
377 }
378
379 return -1 ;
380
381}
382
383//____________________________________________________________________________
384Float_t AliEMCALDigit::GetDEParent(Int_t index) const
385{
386 // retrieves the parent particle energy contribution
387 // given its index in the list
388 if ( (index <= fNiparent) && (index > 0)){
389 return fDEParent[index-1] ;
390 }
391
392 return 0;
393}
394
395//____________________________________________________________________________
396void AliEMCALDigit::ShiftPrimary(Int_t shift){
397 //shifts primary number to BIG offset, to separate primary in different TreeK
398 Int_t index ;
399 for(index = 0; index <fNprimary; index++ ){
400 fPrimary[index] = fPrimary[index]+ shift * 10000000 ;}
401 for(index =0; index <fNiparent; index++){
402 fIparent[index] = fIparent[index] + shift * 10000000 ;}
403}
404//____________________________________________________________________________
405Bool_t AliEMCALDigit::operator==(AliEMCALDigit const & digit) const
406{
407 // Two digits are equal if they have the same Id
408
409 if ( fId == digit.fId )
410 return kTRUE ;
411 else
412 return kFALSE ;
413}
414
415//____________________________________________________________________________
416AliEMCALDigit AliEMCALDigit::operator+(const AliEMCALDigit &digit)
417{
418 // Adds the amplitude of digits and completes the list of primary particles
419 // if amplitude is larger than
420
421 fAmpFloat += digit.fAmpFloat ;
422 for (Int_t i=0; i < fNSamples ; i++) fSamples[i] += digit.fSamples[i];
423 for (Int_t i=0; i < fNSamplesHG; i++) fSamplesHG[i] += digit.fSamplesHG[i];
424
425 fAmp += digit.fAmp ;
426 if(fTime > digit.fTime)
427 fTime = digit.fTime ;
428 if (digit.fTimeR < fTimeR)
429 fTimeR = digit.fTimeR ;
430
431 Int_t max1 = fNprimary ;
432 Int_t max2 = fNiparent ;
433 Int_t index ;
434 for (index = 0 ; index < digit.fNprimary ; index++){
435 Bool_t newPrim = kTRUE ;
436 Int_t old ;
437 for ( old = 0 ; (old < max1) && newPrim; old++) { //already have this primary?
438 if(fPrimary[old] == digit.fPrimary[index]) {
439 newPrim = kFALSE;
440 fDEPrimary[old] += digit.fDEPrimary[index];
441 }
442 }
443 if (newPrim) {
444 if(max1<fNMaxPrimary){
445 fPrimary[max1] = digit.fPrimary[index] ;
446 fDEPrimary[max1] = digit.fDEPrimary[index] ;
447 fNprimary++ ;
448 max1++;
449 }
450 if(fNprimary==fNMaxPrimary) {
451
452 TString mess = " NMaxPrimary = " ;
453 mess += fNMaxPrimary ;
454 mess += " is too small" ;
455 AliFatal(mess.Data()) ;
456
457 }
458 }
459 }
460
461 for (index = 0 ; index < digit.fNiparent ; index++){
462 Bool_t newParent = kTRUE ;
463 Int_t old ;
464 for ( old = 0 ; (old < max2) && newParent; old++) { //already have this primary?
465 if(fIparent[old] == digit.fIparent[index]) {
466 newParent = kFALSE;
467 fDEParent[old] += digit.fDEParent[index];
468 }
469 }
470 if(newParent){
471 if(max2<fNMaxiparent) {
472 fIparent[max2] = digit.fIparent[index] ;
473 fDEParent[max2] = digit.fDEParent[index] ;
474 fNiparent++ ;
475 max2++;
476 }
477 if(fNiparent==fNMaxiparent) {
478
479 TString mess = " NMaxiparent = " ;
480 mess += fNMaxiparent ;
481 mess += " is too small" ;
482 AliFatal(mess.Data()) ;
483
484 }
485 }
486 }
487
488 return *this ;
489}
490
491//____________________________________________________________________________
492AliEMCALDigit AliEMCALDigit::operator*(Float_t factor)
493{
494 // Multiplies the amplitude by a factor
495
496 //Float_t tempo = static_cast<Float_t>(fAmp) ;
497 //tempo *= factor ;
498 //fAmp = static_cast<Int_t>(TMath::Floor(tempo)) ;
499
500 fAmpFloat *= factor;
501 for (Int_t i=0; i < fNSamples ; i++) fSamples[i] = Int_t(factor*fSamples[i]);
502 for (Int_t i=0; i < fNSamplesHG; i++) fSamplesHG[i] = Int_t(factor*fSamplesHG[i]);
503
504 for(Int_t i=0; i < fNprimary; i++)
505 fDEPrimary[i] *= factor;
506 for(Int_t i=0; i < fNiparent; i++)
507 fDEParent[i] *= factor;
508
509 return *this ;
510}
511
512//____________________________________________________________________________
513ostream& operator << ( ostream& out , const AliEMCALDigit & digit)
514{
515 // Prints the data of the digit
516
517 out << "ID " << digit.fId << " Energy = " << digit.fAmp << " Time = " << digit.fTime << endl ;
518 for(Int_t i=0;i<digit.fNprimary;i++)
519 out << "Primary " << i+1 << " = " << digit.fPrimary[i]
520 << " : DE " << digit.fDEPrimary[i] << endl ;
521
522 for(Int_t j=0;j<digit.fNiparent;j++)
523 out << "Iparent " << j+1 << " = " << digit.fIparent[j]
524 << " : DE " << digit.fDEParent[j] << endl ;
525 out << "Position in list = " << digit.fIndexInList << endl ;
526 return out ;
527}
528
529//____________________________________________________________________________
530void AliEMCALDigit::Print(const Option_t* /*opt*/) const
531{
532 //Print
533 printf("===\nDigit id: %4d / Energy %2.3f ; Time %e ; Time samples %d ; Chi2 %2.3f, NDF %d, Type? %d \n",
534 fId, fAmpFloat,fTime, fNSamples, fChi2, fNDF, fDigitType);
535 if(fDigitType==kTrigger){
536 printf("FALTRO: ");
537 for (Int_t i=0; i < GetNFALTROSamples(); i++)
538 {
539 Int_t timeBin, amp;
540 GetFALTROSample(i, timeBin, amp);
541 printf(" (%d,%d) ",timeBin,amp);
542 }
543 printf("\n");
544 }//trigger
545 else{
546 printf("ALTRO, Low Gain: ");
547 for (Int_t i=0; i < GetNALTROSamplesLG(); i++)
548 {
549 Int_t timeBin, amp;
550 GetALTROSampleLG(i, timeBin, amp);
551 printf(" (%d,%d) ",timeBin,amp);
552 }
553 printf("\n");
554 printf("ALTRO, High Gain: ");
555 for (Int_t i=0; i < GetNALTROSamplesHG(); i++)
556 {
557 Int_t timeBin, amp;
558 GetALTROSampleHG(i, timeBin, amp);
559 printf(" (%d,%d) ",timeBin,amp);
560 }
561 printf("\n");
562 }//trigger
563
564}
565
566
567