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