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