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