Adding a reminder for coders
[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   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 //____________________________________________________________________________
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     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 //____________________________________________________________________________
149 AliEMCALDigit::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 //____________________________________________________________________________
218 AliEMCALDigit::~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 //____________________________________________________________________________
233 void 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 //____________________________________________________________________________
246 Int_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 //____________________________________________________________________________
269 Float_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 //____________________________________________________________________________
281 Float_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 //____________________________________________________________________________
293 Bool_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 //____________________________________________________________________________
304 Bool_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 //____________________________________________________________________________
316 void 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 //____________________________________________________________________________
325 void 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 //____________________________________________________________________________
335 Bool_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 //____________________________________________________________________________
348 Int_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 //____________________________________________________________________________
359 Float_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 //____________________________________________________________________________
372 Int_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 //____________________________________________________________________________
384 Float_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 //____________________________________________________________________________
396 void 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 //____________________________________________________________________________
405 Bool_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 //____________________________________________________________________________
416 AliEMCALDigit 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 //____________________________________________________________________________
492 AliEMCALDigit 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 //____________________________________________________________________________
513 ostream& 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 //____________________________________________________________________________
530 void 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