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