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