several changes: 1) check for overflow before fitting, not using fit results 2) keep...
[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   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         
183   fPrimary = new Int_t[fNMaxPrimary] ;  
184   fDEPrimary = new Float_t[fNMaxPrimary] ;
185   fIparent = new Int_t[fNMaxiparent] ;
186   fDEParent = new Float_t[fNMaxiparent] ;
187   Int_t i ;
188   for ( i = 0; i < fNMaxPrimary ; i++) {
189     fPrimary[i]  = digit.fPrimary[i] ;
190     fDEPrimary[i]  = digit.fDEPrimary[i] ;
191   }
192   Int_t j ;
193   for (j = 0; j< fNMaxiparent ; j++) {
194     fIparent[j]  = digit.fIparent[j] ;
195     fDEParent[j]  = digit.fDEParent[j] ;
196   }
197 }
198
199 //____________________________________________________________________________
200 AliEMCALDigit::~AliEMCALDigit() 
201 {
202   // Delete array of primiries if any
203     delete [] fPrimary ;
204     delete [] fDEPrimary ;
205     delete [] fIparent ; 
206     delete [] fDEParent ;
207         delete [] fSamples;
208         delete [] fSamplesHG;
209
210 }
211
212 //____________________________________________________________________________
213 Int_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 //____________________________________________________________________________
236 Float_t AliEMCALDigit::GetEta() const
237
238   //return pseudorapidity for this digit
239   // should be change in EMCALGeometry - 19-nov-04
240   Float_t eta=-10., phi=-10.;
241   Int_t id = GetId();
242   const AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
243   g->EtaPhiFromIndex(id,eta,phi);
244   return eta ;
245 }
246
247 //____________________________________________________________________________
248 Float_t AliEMCALDigit::GetPhi() const
249
250   //return phi coordinate of digit
251   // should be change in EMCALGeometry - 19-nov-04
252   Float_t eta=-10., phi=-10.;
253   Int_t id = GetId();
254   const AliEMCALGeometry *g = AliEMCALGeometry::GetInstance();
255   g->EtaPhiFromIndex(id,eta,phi);
256   return phi ;
257 }
258
259 //____________________________________________________________________________
260 Bool_t AliEMCALDigit::GetFALTROSample(const Int_t iSample, Int_t& timeBin, Int_t& amp) const
261 {
262         //Get FALTRO sample in time bin iSample
263         if (iSample >= fNSamples || iSample < 0 || fDigitType==kTrigger) return kFALSE;
264         
265         amp     =  fSamples[iSample] & 0xFFF;
266         timeBin = (fSamples[iSample] >> 12) & 0xFF;
267         
268         return kTRUE;
269 }
270 //____________________________________________________________________________
271 Bool_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
274         if (iSample >= fNSamples || iSample < 0 || fDigitType==kLG) return kFALSE;
275         
276         amp     =  fSamples[iSample] & 0xFFF;
277         timeBin = (fSamples[iSample] >> 12) & 0xFF;
278         
279         return kTRUE;
280 }
281
282 //____________________________________________________________________________
283 void 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 //____________________________________________________________________________
292 void 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 //____________________________________________________________________________
302 Bool_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
305         if (iSample >= fNSamplesHG || iSample < 0 || fDigitType==kHG) return kFALSE;
306         
307         amp     =  fSamplesHG[iSample] & 0xFFF;
308         timeBin = (fSamplesHG[iSample] >> 12) & 0xFF;
309         
310         return kTRUE;
311 }
312
313
314 //____________________________________________________________________________
315 Int_t AliEMCALDigit::GetPrimary(Int_t index) const
316 {
317   // retrieves the primary particle number given its index in the list 
318   if ( (index <= fNprimary) && (index > 0)){
319     return fPrimary[index-1] ;
320   } 
321
322   return -1 ; 
323 }
324
325 //____________________________________________________________________________
326 Float_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 ; 
335   
336 }
337
338 //____________________________________________________________________________
339 Int_t AliEMCALDigit::GetIparent(Int_t index) const
340 {
341   // retrieves the primary particle number given its index in the list 
342   if ( index <= fNiparent && index > 0){
343     return fIparent[index-1] ;
344   } 
345
346   return -1 ; 
347   
348 }
349
350 //____________________________________________________________________________
351 Float_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 //____________________________________________________________________________
363 void AliEMCALDigit::ShiftPrimary(Int_t shift){
364   //shifts primary number to BIG offset, to separate primary in different TreeK
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 //____________________________________________________________________________
372 Bool_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 //____________________________________________________________________________
383 AliEMCALDigit AliEMCALDigit::operator+(const AliEMCALDigit &digit) 
384 {
385   // Adds the amplitude of digits and completes the list of primary particles
386   // if amplitude is larger than 
387   
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 ;
393   if(fTime > digit.fTime)
394     fTime = digit.fTime ;
395   if (digit.fTimeR < fTimeR)
396     fTimeR = digit.fTimeR ; 
397
398   Int_t max1 = fNprimary ; 
399   Int_t max2 = fNiparent ;  
400   Int_t index ; 
401   for (index = 0 ; index < digit.fNprimary  ; index++){
402     Bool_t newPrim = kTRUE ;
403     Int_t old ;
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       }
409     }
410     if (newPrim) {
411       if(max1<fNMaxPrimary){ 
412         fPrimary[max1] = digit.fPrimary[index] ; 
413         fDEPrimary[max1] = digit.fDEPrimary[index] ; 
414         fNprimary++ ;
415         max1++;
416       }
417       if(fNprimary==fNMaxPrimary) {
418         
419         TString mess = " NMaxPrimary  =  " ; 
420         mess += fNMaxPrimary ; 
421         mess += " is too small" ; 
422         Fatal("AliEMCALDigit::Operator+ -->" , mess.Data()) ; 
423
424       }
425     }
426   }
427   
428   for (index = 0 ; index < digit.fNiparent ; index++){
429     Bool_t newParent = kTRUE ;
430     Int_t old ;
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       }
436     }
437     if(newParent){
438       if(max2<fNMaxiparent) { 
439         fIparent[max2] = digit.fIparent[index] ; 
440         fDEParent[max2] = digit.fDEParent[index] ; 
441         fNiparent++ ;
442         max2++;
443       }
444       if(fNiparent==fNMaxiparent) {
445         
446         TString mess = " NMaxiparent  =  " ; 
447         mess += fNMaxiparent ; 
448         mess += " is too small" ; 
449         Fatal("AliEMCALDigit::Operator+ -->", mess.Data()) ; 
450
451       }
452     }
453   }
454   
455   return *this ;
456 }
457
458 //____________________________________________________________________________
459 AliEMCALDigit AliEMCALDigit::operator*(Float_t factor) 
460 {
461   // Multiplies the amplitude by a factor
462   
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
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
476   return *this ;
477 }
478
479 //____________________________________________________________________________
480 ostream& operator << ( ostream& out , const AliEMCALDigit & digit)
481 {
482   // Prints the data of the digit
483   
484   out << "ID " << digit.fId << " Energy = " << digit.fAmp <<  " Time = " << digit.fTime << endl ; 
485   Int_t i,j ;
486   for(i=0;i<digit.fNprimary;i++) 
487     out << "Primary " << i+1 << " = " << digit.fPrimary[i] 
488         << " : DE " << digit.fDEPrimary[i] << endl ;
489    
490   for(j=0;j<digit.fNiparent;j++)
491     out << "Iparent " << j+1 << " = " << digit.fIparent[j] 
492         << " : DE " << digit.fDEParent[j] << endl ;
493   out << "Position in list = " << digit.fIndexInList << endl ; 
494   return out ;
495 }
496
497 //____________________________________________________________________________
498 void AliEMCALDigit::Print(const Option_t* /*opt*/) const
499 {
500     //Print
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){
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
535