PHOS/AliPHOSEmcRecPoint.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 //_________________________________________________________________________
  17 // Rec Point in the PHOS EM calorimeter
  18 //*-- Author : Dmitri Peressounko RRC KI
  19 //////////////////////////////////////////////////////////////////////////////
  20
  21 // --- ROOT system ---
  22 #include "TPad.h"
  23 #include "TH2.h"
  24 #include "TMath.h"
  25 #include "TCanvas.h"
  26
  27 // --- Standard library ---
  28
  29 #include <iostream>
  30
  31 // --- AliRoot header files ---
  32
  33 #include "AliPHOSGeometry.h"
  34 #include "AliPHOSEmcRecPoint.h"
  35 #include "AliRun.h"
  36
  37 ClassImp(AliPHOSEmcRecPoint)
  38
  39 //____________________________________________________________________________
  40 AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
  41   : AliPHOSRecPoint()
  42 {
  43   // ctor
  44
  45   fMulDigit   = 0 ;
  46   fAmp   = 0. ;
  47   fEnergyList = new Float_t[fMaxDigit];
  48   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
  49   fDelta     =  PHOSGeom->GetCrystalSize(0) ;
  50   fW0        = W0 ;
  51   fLocMaxCut = LocMaxCut ;
  52   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
  53 }
  54
  55 //____________________________________________________________________________
  56 AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint()
  57 {
  58   // dtor
  59 }
  60
  61 //____________________________________________________________________________
  62 void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
  63 {
  64   // adds a digit to the digits list
  65   // and accumulates the total amplitude and the multiplicity
  66
  67   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
  68     fMaxDigit*=2 ;
  69     int * tempo = new ( int[fMaxDigit] ) ;
  70     Float_t * tempoE =  new ( Float_t[fMaxDigit] ) ;
  71
  72     Int_t index ;
  73     for ( index = 0 ; index < fMulDigit ; index++ ){
  74       tempo[index] = fDigitsList[index] ;
  75       tempoE[index] = fEnergyList[index] ;
  76     }
  77
  78     delete [] fDigitsList ;
  79     fDigitsList =  new ( int[fMaxDigit] ) ;
  80
  81     delete [] fEnergyList ;
  82     fEnergyList =  new ( Float_t[fMaxDigit] ) ;
  83
  84     for ( index = 0 ; index < fMulDigit ; index++ ){
  85       fDigitsList[index] = tempo[index] ;
  86       fEnergyList[index] = tempoE[index] ;
  87     }
  88
  89     delete [] tempo ;
  90     delete [] tempoE ;
  91   } // if
  92
  93   fDigitsList[fMulDigit]   =  (int) &digit  ;
  94   fEnergyList[fMulDigit++] = Energy ;
  95   fAmp += Energy ;
  96 }
  97
  98 //____________________________________________________________________________
  99 Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 )
 100 {
 101
 102   Bool_t aren = kFALSE ;
 103
 104   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
 105   Int_t relid1[4] ;
 106   PHOSGeom->AbsToRelNumbering(digit1->GetId(), relid1) ;
 107
 108   Int_t relid2[4] ;
 109   PHOSGeom->AbsToRelNumbering(digit2->GetId(), relid2) ;
 110
 111   Int_t RowDiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 112   Int_t ColDiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 113
 114   if (( ColDiff<=1 )  && ( RowDiff <= 1 ) && (ColDiff+RowDiff > 0))
 115     aren = kTRUE ;
 116
 117   return aren ;
 118 }
 119
 120 //____________________________________________________________________________
 121 Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
 122 {
 123   Int_t rv ;
 124
 125   AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ;
 126
 127
 128   Int_t PHOSMod1 = this->GetPHOSMod() ;
 129   Int_t PHOSMod2 = clu->GetPHOSMod() ;
 130
 131   TVector3 LocPos1;
 132   this->GetLocalPosition(LocPos1) ;
 133   TVector3 LocPos2;
 134   clu->GetLocalPosition(LocPos2) ;
 135
 136   if(PHOSMod1 == PHOSMod2 ) {
 137     Int_t rowdif = (Int_t)TMath::Ceil(LocPos1.X()/fDelta)-(Int_t)TMath::Ceil(LocPos2.X()/fDelta) ;
 138     if (rowdif> 0)
 139       rv = -1 ;
 140      else if(rowdif < 0)
 141        rv = 1 ;
 142     else if(LocPos1.Z()>LocPos2.Z())
 143       rv = -1 ;
 144     else
 145       rv = 1 ;
 146      }
 147
 148   else {
 149     if(PHOSMod1 < PHOSMod2 )
 150       rv = -1 ;
 151     else
 152       rv = 1 ;
 153   }
 154
 155   return rv ;
 156 }
 157
 158 //______________________________________________________________________________
 159 void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
 160 {
 161   //Execute action corresponding to one event
 162   //  This member function is called when a AliPHOSRecPoint is clicked with the locator
 163   //
 164   //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on
 165   //  and switched off when the mouse button is released.
 166   //
 167
 168   //   static Int_t pxold, pyold;
 169
 170    static TGraph *  DigitGraph = 0 ;
 171
 172    if (!gPad->IsEditable()) return;
 173
 174    TH2F * Histo = 0 ;
 175    TCanvas * HistoCanvas ;
 176
 177    switch (event) {
 178
 179    case kButton1Down: {
 180      AliPHOSDigit * digit ;
 181      AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
 182      Int_t iDigit;
 183      Int_t relid[4] ;
 184      Float_t xi[fMulDigit] ;
 185      Float_t zi[fMulDigit] ;
 186
 187      // create the histogram for the single cluster
 188      // 1. gets histogram boundaries
 189      Float_t ximax = -999. ;
 190      Float_t zimax = -999. ;
 191      Float_t ximin = 999. ;
 192      Float_t zimin = 999. ;
 193
 194      for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 195        digit = (AliPHOSDigit *) fDigitsList[iDigit];
 196        PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 197        PHOSGeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
 198        if ( xi[iDigit] > ximax )
 199          ximax = xi[iDigit] ;
 200        if ( xi[iDigit] < ximin )
 201          ximin = xi[iDigit] ;
 202        if ( zi[iDigit] > zimax )
 203          zimax = zi[iDigit] ;
 204        if ( zi[iDigit] < zimin )
 205          zimin = zi[iDigit] ;
 206      }
 207      ximax += PHOSGeom->GetCrystalSize(0) / 2. ;
 208      zimax += PHOSGeom->GetCrystalSize(2) / 2. ;
 209      ximin -= PHOSGeom->GetCrystalSize(0) / 2. ;
 210      zimin -= PHOSGeom->GetCrystalSize(2) / 2. ;
 211      Int_t xdim = (int)( (ximax - ximin ) / PHOSGeom->GetCrystalSize(0) + 0.5  ) ;
 212      Int_t zdim = (int)( (zimax - zimin ) / PHOSGeom->GetCrystalSize(2) + 0.5 ) ;
 213
 214      // 2. gets the histogram title
 215
 216      Text_t title[100] ;
 217      sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
 218
 219      if (!Histo) {
 220        delete Histo ;
 221        Histo = 0 ;
 222      }
 223      Histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
 224
 225      Float_t x, z ;
 226      for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 227        digit = (AliPHOSDigit *) fDigitsList[iDigit];
 228        PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 229        PHOSGeom->RelPosInModule(relid, x, z);
 230        Histo->Fill(x, z, fEnergyList[iDigit] ) ;
 231      }
 232
 233      if (!DigitGraph) {
 234        DigitGraph = new TGraph(fMulDigit,xi,zi);
 235        DigitGraph-> SetMarkerStyle(5) ;
 236        DigitGraph-> SetMarkerSize(1.) ;
 237        DigitGraph-> SetMarkerColor(1) ;
 238        DigitGraph-> Paint("P") ;
 239      }
 240
 241      Print() ;
 242      HistoCanvas = new TCanvas("cluser", "a single cluster", 600, 500) ;
 243      HistoCanvas->Draw() ;
 244      Histo->Draw("lego1") ;
 245
 246      break;
 247    }
 248
 249    case kButton1Up:
 250      if (DigitGraph) {
 251        delete DigitGraph  ;
 252        DigitGraph = 0 ;
 253      }
 254      break;
 255
 256    }
 257 }
 258
 259 //____________________________________________________________________________
 260 Float_t  AliPHOSEmcRecPoint::GetDispersion()
 261 {
 262   Float_t D    = 0 ;
 263   Float_t wtot = 0 ;
 264
 265   TVector3 LocPos;
 266   GetLocalPosition(LocPos);
 267   Float_t x = LocPos.X() ;
 268   Float_t z = LocPos.Z() ;
 269   //  Int_t i = GetPHOSMod() ;
 270
 271   AliPHOSDigit * digit ;
 272   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
 273
 274   Int_t iDigit;
 275   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 276     digit = (AliPHOSDigit *) fDigitsList[iDigit];
 277     Int_t relid[4] ;
 278     Float_t xi ;
 279     Float_t zi ;
 280     PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 281     PHOSGeom->RelPosInModule(relid, xi, zi);
 282     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 283     D += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ;
 284     wtot+=w ;
 285   }
 286
 287   D /= wtot ;
 288
 289   return TMath::Sqrt(D) ;
 290 }
 291
 292 //____________________________________________________________________________
 293 void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
 294 {
 295   Float_t wtot = 0. ;
 296   Float_t x    = 0.;
 297   Float_t z    = 0.;
 298   Float_t Dxx  = 0.;
 299   Float_t Dzz  = 0.;
 300   Float_t Dxz  = 0.;
 301
 302   AliPHOSDigit * digit ;
 303   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
 304   Int_t iDigit;
 305
 306   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 307     digit = (AliPHOSDigit *) fDigitsList[iDigit];
 308     Int_t relid[4] ;
 309     Float_t xi ;
 310     Float_t zi ;
 311     PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 312     PHOSGeom->RelPosInModule(relid, xi, zi);
 313     Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 314     Dxx  += w * xi * xi ;
 315     x    += w * xi ;
 316     Dzz  += w * zi * zi ;
 317     z    += w * zi ;
 318     Dxz  += w * xi * zi ;
 319     wtot += w ;
 320   }
 321
 322   Dxx /= wtot ;
 323   x   /= wtot ;
 324   Dxx -= x * x ;
 325   Dzz /= wtot ;
 326   z   /= wtot ;
 327   Dzz -= z * z ;
 328   Dxz /= wtot ;
 329   Dxz -= x * z ;
 330
 331   lambda[0] = TMath::Sqrt( 0.5 * (Dxx + Dzz) + TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
 332   lambda[1] = TMath::Sqrt( 0.5 * (Dxx + Dzz) - TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
 333 }
 334
 335 //____________________________________________________________________________
 336 Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
 337 {
 338   Float_t menergy = 0. ;
 339
 340   Int_t iDigit;
 341
 342   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 343
 344     if(fEnergyList[iDigit] > menergy)
 345       menergy = fEnergyList[iDigit] ;
 346   }
 347   return menergy ;
 348 }
 349
 350 //____________________________________________________________________________
 351 Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H)
 352 {
 353   Int_t multipl   = 0 ;
 354   Int_t iDigit ;
 355   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 356
 357     if(fEnergyList[iDigit] > H * fAmp)
 358       multipl++ ;
 359   }
 360   return multipl ;
 361 }
 362
 363 //____________________________________________________________________________
 364 Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy)
 365 {
 366   AliPHOSDigit * digit ;
 367   AliPHOSDigit * digitN ;
 368
 369
 370   Int_t iDigitN ;
 371   Int_t iDigit ;
 372
 373   for(iDigit = 0; iDigit < fMulDigit; iDigit++){
 374     maxAt[iDigit] = fDigitsList[iDigit] ;
 375   }
 376
 377   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
 378     if(maxAt[iDigit] != -1) {
 379       digit = (AliPHOSDigit *) maxAt[iDigit] ;
 380
 381       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
 382         digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ;
 383
 384         if ( AreNeighbours(digit, digitN) ) {
 385           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
 386             maxAt[iDigitN] = -1 ;
 387             // but may be digit too is not local max ?
 388             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut)
 389               maxAt[iDigit] = -1 ;
 390           }
 391           else {
 392             maxAt[iDigit] = -1 ;
 393             // but may be digitN too is not local max ?
 394             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut)
 395               maxAt[iDigitN] = -1 ;
 396           }
 397         } // if Areneighbours
 398       } // while digitN
 399     } // slot not empty
 400   } // while digit
 401
 402   iDigitN = 0 ;
 403   for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
 404     if(maxAt[iDigit] != -1){
 405       maxAt[iDigitN] = maxAt[iDigit] ;
 406       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
 407       iDigitN++ ;
 408     }
 409   }
 410   return iDigitN ;
 411 }
 412
 413 //____________________________________________________________________________
 414 void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
 415 {
 416   if( fLocPos.X() < 1000000.) { // already evaluated
 417    LPos = fLocPos ;
 418    return ;
 419   }
 420
 421   Float_t wtot = 0. ;
 422
 423   Int_t relid[4] ;
 424
 425   Float_t x = 0. ;
 426   Float_t z = 0. ;
 427
 428   AliPHOSDigit * digit ;
 429
 430   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
 431
 432   Int_t iDigit;
 433
 434
 435   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 436     digit = (AliPHOSDigit *) fDigitsList[iDigit];
 437
 438     Float_t xi ;
 439     Float_t zi ;
 440     PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 441     PHOSGeom->RelPosInModule(relid, xi, zi);
 442     Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
 443     x    += xi * w ;
 444     z    += zi * w ;
 445     wtot += w ;
 446
 447   }
 448
 449   x /= wtot ;
 450   z /= wtot ;
 451   fLocPos.SetX(x)  ;
 452   fLocPos.SetY(0.) ;
 453   fLocPos.SetZ(z)  ;
 454
 455   LPos = fLocPos ;
 456 }
 457
 458 // //____________________________________________________________________________
 459 // AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu)
 460 // {
 461 //   int * DL = Clu.GetDigitsList() ;
 462
 463 //  if(fDigitsList)
 464 //     delete fDigitsList  ;
 465
 466 //   AliPHOSDigit * digit ;
 467
 468 //   Int_t iDigit;
 469
 470 //   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 471 //     digit = (AliPHOSDigit *) DL[iDigit];
 472 //     AddDigit(*digit) ;
 473 //   }
 474
 475 //   fAmp       = Clu.GetTotalEnergy() ;
 476 //   fGeom      = Clu.GetGeom() ;
 477 //   TVector3 LocPos;
 478 //   Clu.GetLocalPosition(LocPos) ;
 479 //   fLocPos    = LocPos;
 480 //   fMulDigit       = Clu.GetMultiplicity() ;
 481 //   fMaxDigit  = Clu.GetMaximumMultiplicity() ;
 482 //   fPHOSMod   = Clu.GetPHOSMod() ;
 483 //   fW0        = Clu.GetLogWeightCut() ;
 484 //   fDelta     = Clu.GetDelta() ;
 485 //   fLocMaxCut = Clu.GetLocMaxCut() ;
 486
 487 //   delete DL ;
 488
 489 //   return *this ;
 490 // }
 491
 492 //____________________________________________________________________________
 493 void AliPHOSEmcRecPoint::Print(Option_t * option)
 494 {
 495   cout << "AliPHOSEmcRecPoint: " << endl ;
 496
 497   AliPHOSDigit * digit ;
 498   Int_t iDigit;
 499   AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
 500
 501   Float_t xi ;
 502   Float_t zi ;
 503   Int_t relid[4] ;
 504
 505   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 506     digit = (AliPHOSDigit *) fDigitsList[iDigit];
 507     PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
 508     PHOSGeom->RelPosInModule(relid, xi, zi);
 509     cout << " Id = " << digit->GetId() ;
 510     cout << "   module  = " << relid[0] ;
 511     cout << "   x  = " << xi ;
 512     cout << "   z  = " << zi ;
 513     cout << "   Energy = " << fEnergyList[iDigit] << endl ;
 514   }
 515   cout << "       Multiplicity    = " << fMulDigit  << endl ;
 516   cout << "       Cluster Energy  = " << fAmp << endl ;
 517
 518 }
 519
 520 //______________________________________________________________________________
 521 void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
 522 {
 523    // Stream an object of class AliPHOSEmcRecPoint.
 524
 525    if (R__b.IsReading()) {
 526       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 527       AliPHOSRecPoint::Streamer(R__b);
 528       R__b >> fDelta;
 529       R__b >> fLocMaxCut;
 530       R__b.ReadArray(fEnergyList);
 531       R__b >> fW0;
 532    } else {
 533       R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
 534       AliPHOSRecPoint::Streamer(R__b);
 535       R__b << fDelta;
 536       R__b << fLocMaxCut;
 537       R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
 538       R__b << fW0;
 539    }
 540 }
 541
 542