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   fTriggerDigit(kFALSE)
  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, Bool_t trigger, Int_t index, 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(0.),
 104     fNDF(0),
 105         fTriggerDigit(trigger)
 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         fTriggerDigit(digit.fTriggerDigit)
 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 || !fTriggerDigit) 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 || fTriggerDigit) 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 || fTriggerDigit) 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, Trigger? %d \n",
 502                    fId, fAmpFloat,fTime, fNSamples, fChi2, fNDF, fTriggerDigit);
 503         if(fTriggerDigit){
 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