PHOS/PHOSbase/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 /* $Id$ */
  17
  18 /* History of cvs commits:
  19  *
  20  * $Log$
  21  * Revision 1.59  2007/10/18 15:12:22  kharlov
  22  * Moved MakePrimary to EMCRecPoint to rpduce correct order of primaries
  23  *
  24  * Revision 1.58  2007/04/16 09:03:37  kharlov
  25  * Incedent angle correction fixed
  26  *
  27  * Revision 1.57  2007/04/05 10:18:58  policheh
  28  * Introduced distance to nearest bad crystal.
  29  *
  30  * Revision 1.56  2007/03/06 06:47:28  kharlov
  31  * DP:Possibility to use actual vertex position added
  32  *
  33  * Revision 1.55  2007/01/19 20:31:19  kharlov
  34  * Improved formatting for Print()
  35  */
  36
  37 //_________________________________________________________________________
  38 //  RecPoint implementation for PHOS-EMC
  39 //  An EmcRecPoint is a cluster of digits
  40 //--
  41 //-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
  42
  43
  44 // --- ROOT system ---
  45 #include "TH2.h"
  46 #include "TMath.h"
  47 #include "TCanvas.h"
  48 #include "TGraph.h"
  49
  50 // --- Standard library ---
  51
  52 // --- AliRoot header files ---
  53 #include "AliLog.h"
  54 #include "AliPHOSLoader.h"
  55 #include "AliGenerator.h"
  56 #include "AliPHOSGeometry.h"
  57 #include "AliPHOSDigit.h"
  58 #include "AliPHOSEmcRecPoint.h"
  59 #include "AliPHOSReconstructor.h"
  60
  61 ClassImp(AliPHOSEmcRecPoint)
  62
  63 Long64_t AliPHOSEmcRecPoint::fgInstCount=0;
  64
  65
  66 //____________________________________________________________________________
  67 AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() :
  68   AliPHOSRecPoint(),
  69   fCoreEnergy(0.), fDispersion(0.),
  70   fEnergyList(0), fTime(-1.), fNExMax(0),
  71   fM2x(0.), fM2z(0.), fM3x(0.), fM4z(0.),
  72   fPhixe(0.), fDistToBadCrystal(-1),fDebug(0)
  73   ,fInstCount(0)
  74 {
  75   // ctor
  76   fMulDigit   = 0 ;
  77   fAmp   = 0. ;
  78   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
  79
  80   fLambda[0] = 0.;
  81   fLambda[1] = 0.;
  82   fInstCount=fgInstCount++;
  83   if (gDebug==-10) AliInfo(Form("0Create instance %lld",fInstCount));
  84 }
  85
  86 //____________________________________________________________________________
  87 AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const char * opt) :
  88   AliPHOSRecPoint(opt),
  89   fCoreEnergy(0.), fDispersion(0.),
  90   fEnergyList(0), fTime(-1.), fNExMax(0),
  91   fM2x(0.), fM2z(0.), fM3x(0.), fM4z(0.),
  92   fPhixe(0.), fDistToBadCrystal(-1), fDebug(0)
  93   ,fInstCount(0)
  94 {
  95   // ctor
  96   fMulDigit   = 0 ;
  97   fAmp   = 0. ;
  98   fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
  99
 100   fLambda[0] = 0.;
 101   fLambda[1] = 0.;
 102   fInstCount=fgInstCount++;
 103    if (gDebug==-10) AliInfo(Form("1Create instance %lld",fInstCount));
 104 }
 105
 106 //____________________________________________________________________________
 107 AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const AliPHOSEmcRecPoint & rp) :
 108   AliPHOSRecPoint(rp),
 109   fCoreEnergy(rp.fCoreEnergy), fDispersion(rp.fDispersion),
 110   fEnergyList(0), fTime(rp.fTime), fNExMax(rp.fNExMax),
 111   fM2x(rp.fM2x), fM2z(rp.fM2z), fM3x(rp.fM3x), fM4z(rp.fM4z),
 112   fPhixe(rp.fPhixe), fDistToBadCrystal(rp.fDistToBadCrystal), fDebug(rp.fDebug)
 113   ,fInstCount(0)
 114 {
 115   // cpy ctor
 116   fMulDigit   = rp.fMulDigit ;
 117   fAmp        = rp.fAmp ;
 118   if (rp.fMulDigit>0) fEnergyList = new Float_t[rp.fMulDigit] ;
 119   for(Int_t index = 0 ; index < fMulDigit ; index++)
 120     fEnergyList[index] = rp.fEnergyList[index] ;
 121
 122   for(Int_t i=0; i<2; i++) {
 123     fLambda[i] = rp.fLambda[i];
 124   }
 125   fInstCount=fgInstCount++;
 126    if (gDebug==-10) AliInfo(Form("2Create instance %lld",fInstCount));
 127 }
 128
 129 //____________________________________________________________________________
 130 AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint()
 131 {
 132   // dtor
 133   if ( fEnergyList )
 134     delete[] fEnergyList ;
 135    if (gDebug==-10) AliInfo(Form("Delete instance %lld (%lld)",fInstCount, fgInstCount));
 136   if (fInstCount>=fgInstCount-1) fgInstCount--;
 137
 138 }
 139
 140 //____________________________________________________________________________
 141 void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy, Float_t time)
 142 {
 143   // Adds a digit to the RecPoint
 144   // and accumulates the total amplitude and the multiplicity
 145
 146   if(fEnergyList == 0)
 147     fEnergyList =  new Float_t[fMaxDigit];
 148
 149   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
 150     fMaxDigit*=2 ;
 151     Int_t * tempo = new Int_t[fMaxDigit];
 152     Float_t * tempoE =  new Float_t[fMaxDigit];
 153
 154     Int_t index ;
 155     for ( index = 0 ; index < fMulDigit ; index++ ){
 156       tempo[index]  = fDigitsList[index] ;
 157       tempoE[index] = fEnergyList[index] ;
 158     }
 159
 160     delete [] fDigitsList ;
 161     fDigitsList =  new Int_t[fMaxDigit];
 162
 163     delete [] fEnergyList ;
 164     fEnergyList =  new Float_t[fMaxDigit];
 165
 166     for ( index = 0 ; index < fMulDigit ; index++ ){
 167       fDigitsList[index] = tempo[index] ;
 168       fEnergyList[index] = tempoE[index] ;
 169     }
 170
 171     delete [] tempo ;
 172     delete [] tempoE ;
 173   } // if
 174
 175   //time
 176   Bool_t isMax=kTRUE ;
 177   for(Int_t index = 0 ; index < fMulDigit ; index++ ){
 178     if(fEnergyList[index]>Energy){
 179       isMax=kFALSE ;
 180       break ;
 181     }
 182   }
 183   if(isMax){
 184     fTime=time ;
 185   }
 186   //Alternative time calculation - still to be validated
 187   // fTime = (fTime*fAmp + time*Energy)/(fAmp+Energy) ;
 188
 189   fDigitsList[fMulDigit]   = digit.GetIndexInList()  ;
 190   fEnergyList[fMulDigit]   = Energy ;
 191   fMulDigit++ ;
 192   fAmp += Energy ;
 193   EvalPHOSMod(&digit) ;
 194 }
 195
 196 //____________________________________________________________________________
 197 Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) const
 198 {
 199   // Tells if (true) or not (false) two digits are neighbors
 200
 201   Bool_t aren = kFALSE ;
 202
 203   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 204
 205   Int_t relid1[4] ;
 206   phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ;
 207
 208   Int_t relid2[4] ;
 209   phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ;
 210
 211   Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
 212   Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
 213
 214   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
 215     aren = kTRUE ;
 216
 217   return aren ;
 218 }
 219
 220 //____________________________________________________________________________
 221 Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const
 222 {
 223   // Compares two RecPoints according to their position in the PHOS modules
 224
 225   const Float_t delta = 1 ; //Width of "Sorting row". If you changibg this
 226                       //value (what is senseless) change as vell delta in
 227                       //AliPHOSTrackSegmentMakerv* and other RecPoints...
 228   Int_t rv ;
 229
 230   AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ;
 231
 232
 233   Int_t phosmod1 = GetPHOSMod() ;
 234   Int_t phosmod2 = clu->GetPHOSMod() ;
 235
 236   TVector3 locpos1;
 237   GetLocalPosition(locpos1) ;
 238   TVector3 locpos2;
 239   clu->GetLocalPosition(locpos2) ;
 240
 241   if(phosmod1 == phosmod2 ) {
 242     Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ;
 243     if (rowdif> 0)
 244       rv = 1 ;
 245      else if(rowdif < 0)
 246        rv = -1 ;
 247     else if(locpos1.Z()>locpos2.Z())
 248       rv = -1 ;
 249     else
 250       rv = 1 ;
 251      }
 252
 253   else {
 254     if(phosmod1 < phosmod2 )
 255       rv = -1 ;
 256     else
 257       rv = 1 ;
 258   }
 259
 260   return rv ;
 261 }
 262 //______________________________________________________________________________
 263 void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) /*const*/
 264 {
 265
 266   // Execute action corresponding to one event
 267   //  This member function is called when a AliPHOSRecPoint is clicked with the locator
 268   //
 269   //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on
 270   //  and switched off when the mouse button is released.
 271
 272
 273   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance();
 274
 275   static TGraph *  digitgraph = 0 ;
 276
 277   if (!gPad->IsEditable()) return;
 278
 279   TH2F * histo = 0 ;
 280   TCanvas * histocanvas ;
 281
 282
 283   //try to get run loader from default event folder
 284   AliRunLoader* rn = AliRunLoader::GetRunLoader(AliConfig::GetDefaultEventFolderName());
 285   if (rn == 0x0)
 286     {
 287       AliError(Form("Cannot find Run Loader in Default Event Folder"));
 288       return;
 289     }
 290   AliPHOSLoader* phosLoader = dynamic_cast<AliPHOSLoader*>(rn->GetLoader("PHOSLoader"));
 291   if (phosLoader == 0x0)
 292     {
 293       AliError(Form("Cannot find PHOS Loader from Run Loader"));
 294       return;
 295     }
 296
 297
 298   const TClonesArray * digits = phosLoader->Digits() ;
 299
 300   switch (event) {
 301
 302   case kButton1Down: {
 303     AliPHOSDigit * digit ;
 304     Int_t iDigit;
 305     Int_t relid[4] ;
 306
 307     const Int_t kMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ;
 308     Float_t * xi = new Float_t[kMulDigit] ;
 309     Float_t * zi = new Float_t[kMulDigit] ;
 310
 311     // create the histogram for the single cluster
 312     // 1. gets histogram boundaries
 313     Float_t ximax = -999. ;
 314     Float_t zimax = -999. ;
 315     Float_t ximin = 999. ;
 316     Float_t zimin = 999. ;
 317
 318     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
 319       digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
 320       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 321       phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
 322       if ( xi[iDigit] > ximax )
 323         ximax = xi[iDigit] ;
 324       if ( xi[iDigit] < ximin )
 325         ximin = xi[iDigit] ;
 326       if ( zi[iDigit] > zimax )
 327         zimax = zi[iDigit] ;
 328       if ( zi[iDigit] < zimin )
 329         zimin = zi[iDigit] ;
 330     }
 331     ximax += phosgeom->GetCrystalSize(0) / 2. ;
 332     zimax += phosgeom->GetCrystalSize(2) / 2. ;
 333     ximin -= phosgeom->GetCrystalSize(0) / 2. ;
 334     zimin -= phosgeom->GetCrystalSize(2) / 2. ;
 335     Int_t xdim = (int)( (ximax - ximin ) / phosgeom->GetCrystalSize(0) + 0.5  ) ;
 336     Int_t zdim = (int)( (zimax - zimin ) / phosgeom->GetCrystalSize(2) + 0.5 ) ;
 337
 338     // 2. gets the histogram title
 339
 340     Text_t title[100] ;
 341     snprintf(title,100,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
 342
 343     if (!histo) {
 344       delete histo ;
 345       histo = 0 ;
 346     }
 347     histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;
 348
 349     Float_t x, z ;
 350     for(iDigit=0; iDigit<kMulDigit; iDigit++) {
 351       digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
 352       phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 353       phosgeom->RelPosInModule(relid, x, z);
 354       histo->Fill(x, z, fEnergyList[iDigit] ) ;
 355     }
 356
 357     if (!digitgraph) {
 358       digitgraph = new TGraph(kMulDigit,xi,zi);
 359       digitgraph-> SetMarkerStyle(5) ;
 360       digitgraph-> SetMarkerSize(1.) ;
 361       digitgraph-> SetMarkerColor(1) ;
 362       digitgraph-> Paint("P") ;
 363     }
 364
 365     //    Print() ;
 366     histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ;
 367     histocanvas->Draw() ;
 368     histo->Draw("lego1") ;
 369
 370     delete[] xi ;
 371     delete[] zi ;
 372
 373     break;
 374   }
 375
 376   case kButton1Up:
 377     if (digitgraph) {
 378       delete digitgraph  ;
 379       digitgraph = 0 ;
 380     }
 381     break;
 382
 383    }
 384 }
 385
 386 //____________________________________________________________________________
 387 void  AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits, TVector3 & /* vInc */)
 388 {
 389   // Calculates the dispersion of the shower at the origine of the RecPoint
 390   //DP: should we correct dispersion for non-perpendicular hit????????
 391
 392   Float_t d    = 0. ;
 393   Float_t wtot = 0. ;
 394
 395   Float_t x = 0.;
 396   Float_t z = 0.;
 397
 398   AliPHOSDigit * digit ;
 399
 400   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance();
 401
 402  // Calculates the center of gravity in the local PHOS-module coordinates
 403
 404   Int_t iDigit;
 405
 406   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 407     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ;
 408     Int_t relid[4] ;
 409     Float_t xi ;
 410     Float_t zi ;
 411     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 412     phosgeom->RelPosInModule(relid, xi, zi);
 413     if (fAmp>0 && fEnergyList[iDigit]>0) {
 414       Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
 415       x    += xi * w ;
 416       z    += zi * w ;
 417       wtot += w ;
 418     }
 419     else
 420       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
 421   }
 422   if (wtot>0) {
 423     x /= wtot ;
 424     z /= wtot ;
 425   }
 426   else
 427     AliError(Form("Wrong weight %f\n", wtot));
 428
 429
 430 // Calculates the dispersion in coordinates
 431   wtot = 0.;
 432   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
 433     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
 434     Int_t relid[4] ;
 435     Float_t xi ;
 436     Float_t zi ;
 437     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 438     phosgeom->RelPosInModule(relid, xi, zi);
 439     if (fAmp>0 && fEnergyList[iDigit]>0) {
 440       Float_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 441       d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ;
 442       wtot+=w ;
 443     }
 444     else
 445       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
 446   }
 447
 448
 449   if (wtot>0) {
 450     d /= wtot ;
 451   }
 452   else
 453     AliError(Form("Wrong weight %f\n", wtot));
 454
 455   fDispersion = 0;
 456   if (d>=0)
 457     fDispersion = TMath::Sqrt(d) ;
 458
 459
 460 }
 461 //______________________________________________________________________________
 462 void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, Float_t coreRadius, TClonesArray * digits)
 463 {
 464   // This function calculates energy in the core,
 465   // i.e. within a radius rad = 3cm around the center. Beyond this radius
 466   // in accordance with shower profile the energy deposition
 467   // should be less than 2%
 468 //DP: non-perpendicular incidence??????????????
 469
 470   Float_t x = 0 ;
 471   Float_t z = 0 ;
 472
 473   AliPHOSDigit * digit ;
 474
 475   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance();
 476
 477   Int_t iDigit;
 478
 479 // Calculates the center of gravity in the local PHOS-module coordinates
 480   Float_t wtot = 0;
 481   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 482     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ;
 483     Int_t relid[4] ;
 484     Float_t xi ;
 485     Float_t zi ;
 486     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 487     phosgeom->RelPosInModule(relid, xi, zi);
 488     if (fAmp>0 && fEnergyList[iDigit]>0) {
 489       Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
 490       x    += xi * w ;
 491       z    += zi * w ;
 492       wtot += w ;
 493     }
 494     else
 495       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
 496   }
 497   if (wtot>0) {
 498     x /= wtot ;
 499     z /= wtot ;
 500   }
 501   else
 502     AliError(Form("Wrong weight %f\n", wtot));
 503
 504
 505   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
 506     digit = (AliPHOSDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
 507     Int_t relid[4] ;
 508     Float_t xi ;
 509     Float_t zi ;
 510     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 511     phosgeom->RelPosInModule(relid, xi, zi);
 512     Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(zi-z)*(zi-z)) ;
 513     if(distance < coreRadius)
 514       fCoreEnergy += fEnergyList[iDigit] ;
 515   }
 516
 517
 518 }
 519
 520 //____________________________________________________________________________
 521 void  AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits, TVector3 & /* vInc */)
 522 {
 523   // Calculates the axis of the shower ellipsoid
 524
 525   Double_t wtot = 0. ;
 526   Double_t x    = 0.;
 527   Double_t z    = 0.;
 528   Double_t dxx  = 0.;
 529   Double_t dzz  = 0.;
 530   Double_t dxz  = 0.;
 531
 532   AliPHOSDigit * digit ;
 533
 534   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance();
 535
 536   Int_t iDigit;
 537
 538
 539   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 540     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
 541     Int_t relid[4] ;
 542     Float_t xi ;
 543     Float_t zi ;
 544     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 545     phosgeom->RelPosInModule(relid, xi, zi);
 546     if (fAmp>0 && fEnergyList[iDigit]>0) {
 547       Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 548       dxx  += w * xi * xi ;
 549       x    += w * xi ;
 550       dzz  += w * zi * zi ;
 551       z    += w * zi ;
 552       dxz  += w * xi * zi ;
 553       wtot += w ;
 554     }
 555     else
 556       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
 557   }
 558   if (wtot>0) {
 559     dxx /= wtot ;
 560     x   /= wtot ;
 561     dxx -= x * x ;
 562     dzz /= wtot ;
 563     z   /= wtot ;
 564     dzz -= z * z ;
 565     dxz /= wtot ;
 566     dxz -= x * z ;
 567
 568 //   //Apply correction due to non-perpendicular incidence
 569 //   Double_t CosX ;
 570 //   Double_t CosZ ;
 571 //   AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
 572 //   Double_t DistanceToIP= (Double_t ) phosgeom->GetIPtoCrystalSurface() ;
 573
 574 //   CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ;
 575 //   CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ;
 576
 577 //   dxx = dxx/(CosX*CosX) ;
 578 //   dzz = dzz/(CosZ*CosZ) ;
 579 //   dxz = dxz/(CosX*CosZ) ;
 580
 581
 582     fLambda[0] =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 583     if(fLambda[0] > 0)
 584       fLambda[0] = TMath::Sqrt(fLambda[0]) ;
 585
 586     fLambda[1] =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 587     if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
 588       fLambda[1] = TMath::Sqrt(fLambda[1]) ;
 589     else
 590       fLambda[1]= 0. ;
 591   }
 592   else {
 593     AliError(Form("Wrong weight %f\n", wtot));
 594     fLambda[0]=fLambda[1]=0.;
 595   }
 596 }
 597
 598 //____________________________________________________________________________
 599 void  AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits, TVector3 & /* vInc */)
 600 {
 601   // Calculate the shower moments in the eigen reference system
 602   // M2x, M2z, M3x, M4z
 603   // Calculate the angle between the shower position vector and the eigen vector
 604
 605   Double_t wtot = 0. ;
 606   Double_t x    = 0.;
 607   Double_t z    = 0.;
 608   Double_t dxx  = 0.;
 609   Double_t dzz  = 0.;
 610   Double_t dxz  = 0.;
 611   Double_t lambda0=0, lambda1=0;
 612
 613   AliPHOSDigit * digit ;
 614
 615   AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
 616
 617   Int_t iDigit;
 618
 619   // 1) Find covariance matrix elements:
 620   //    || dxx dxz ||
 621   //    || dxz dzz ||
 622
 623   Float_t xi ;
 624   Float_t zi ;
 625   Int_t relid[4] ;
 626   Double_t w;
 627   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 628     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
 629     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 630     phosgeom->RelPosInModule(relid, xi, zi);
 631     if (fAmp>0 && fEnergyList[iDigit]>0) {
 632       w     = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 633       x    += w * xi ;
 634       z    += w * zi ;
 635       dxx  += w * xi * xi ;
 636       dzz  += w * zi * zi ;
 637       dxz  += w * xi * zi ;
 638       wtot += w ;
 639     }
 640     else
 641       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
 642
 643   }
 644   if (wtot>0) {
 645     x   /= wtot ;
 646     z   /= wtot ;
 647     dxx /= wtot ;
 648     dzz /= wtot ;
 649     dxz /= wtot ;
 650     dxx -= x * x ;
 651     dzz -= z * z ;
 652     dxz -= x * z ;
 653
 654   // 2) Find covariance matrix eigen values lambda0 and lambda1
 655
 656     lambda0 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 657     lambda1 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 658   }
 659   else {
 660     AliError(Form("Wrong weight %f\n", wtot));
 661     lambda0=lambda1=0.;
 662   }
 663
 664   // 3) Find covariance matrix eigen vectors e0 and e1
 665
 666   TVector2 e0,e1;
 667   if (dxz != 0)
 668     e0.Set(1.,(lambda0-dxx)/dxz);
 669   else
 670     e0.Set(0.,1.);
 671
 672   e0 = e0.Unit();
 673   e1.Set(-e0.Y(),e0.X());
 674
 675   // 4) Rotate cluster tensor from (x,z) to (e0,e1) system
 676   //    and calculate moments M3x and M4z
 677
 678   Float_t cosPhi = e0.X();
 679   Float_t sinPhi = e0.Y();
 680
 681   Float_t xiPHOS ;
 682   Float_t ziPHOS ;
 683   Double_t dx3, dz3, dz4;
 684   wtot = 0.;
 685   x    = 0.;
 686   z    = 0.;
 687   dxx  = 0.;
 688   dzz  = 0.;
 689   dxz  = 0.;
 690   dx3  = 0.;
 691   dz3  = 0.;
 692   dz4  = 0.;
 693   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 694     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit])  ;
 695     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 696     phosgeom->RelPosInModule(relid, xiPHOS, ziPHOS);
 697     xi    = xiPHOS*cosPhi + ziPHOS*sinPhi;
 698     zi    = ziPHOS*cosPhi - xiPHOS*sinPhi;
 699     if (fAmp>0 && fEnergyList[iDigit]>0) {
 700       w     = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 701       x    += w * xi ;
 702       z    += w * zi ;
 703       dxx  += w * xi * xi ;
 704       dzz  += w * zi * zi ;
 705       dxz  += w * xi * zi ;
 706       dx3  += w * xi * xi * xi;
 707       dz3  += w * zi * zi * zi ;
 708       dz4  += w * zi * zi * zi * zi ;
 709       wtot += w ;
 710     }
 711     else
 712       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
 713   }
 714   if (wtot>0) {
 715     x   /= wtot ;
 716     z   /= wtot ;
 717     dxx /= wtot ;
 718     dzz /= wtot ;
 719     dxz /= wtot ;
 720     dx3 /= wtot ;
 721     dz3 /= wtot ;
 722     dz4 /= wtot ;
 723     dx3 += -3*dxx*x + 2*x*x*x;
 724     dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z;
 725     dxx -= x * x ;
 726     dzz -= z * z ;
 727     dxz -= x * z ;
 728   }
 729   else
 730     AliError(Form("Wrong weight %f\n", wtot));
 731
 732   // 5) Find an angle between cluster center vector and eigen vector e0
 733
 734   Float_t phi = 0;
 735   if ( (x*x+z*z) > 0 )
 736     phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z));
 737
 738   fM2x   = lambda0;
 739   fM2z   = lambda1;
 740   fM3x   = dx3;
 741   fM4z   = dz4;
 742   fPhixe = phi;
 743
 744 }
 745 //______________________________________________________________________________
 746 void  AliPHOSEmcRecPoint::EvalPrimaries(TClonesArray * digits)
 747 {
 748   // Constructs the list of primary particles (tracks) which have contributed to this RecPoint
 749
 750   AliPHOSDigit * digit ;
 751   Int_t * tempo    = new Int_t[fMaxTrack] ;
 752
 753   //First find digit with maximal energy deposition and copy its primaries
 754   Float_t emax=0.;
 755   Int_t imaxDigit=0;
 756   for(Int_t id=0; id<GetDigitsMultiplicity(); id++){
 757     if(emax<fEnergyList[id]){
 758       imaxDigit=id ;
 759       emax=fEnergyList[id];
 760     }
 761   }
 762   digit = static_cast<AliPHOSDigit *>(digits->At( fDigitsList[imaxDigit] )) ;
 763   Int_t nprimaries = digit->GetNprimary() ;
 764   if ( nprimaries > fMaxTrack ) {
 765     fMulTrack = - 1 ;
 766     Error("EvalPrimaries", "GetNprimaries ERROR > increase fMaxTrack" ) ;
 767     nprimaries = fMaxTrack; //skip the rest
 768   }
 769   for(fMulTrack=0; fMulTrack<nprimaries ; fMulTrack++){
 770     tempo[fMulTrack] = digit->GetPrimary(fMulTrack+1) ;
 771   }
 772
 773   //Now add other digits contributions
 774   for (Int_t index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
 775     if(index==imaxDigit) //already in
 776       continue ;
 777     digit = static_cast<AliPHOSDigit *>(digits->At( fDigitsList[index] )) ;
 778     nprimaries = digit->GetNprimary() ;
 779     for(Int_t ipr=0; ipr<nprimaries; ipr++){
 780       Int_t iprimary = digit->GetPrimary(ipr+1) ;
 781       Bool_t notIn=1 ;
 782       for(Int_t kndex = 0 ; (kndex < fMulTrack)&& notIn ; kndex++ ) { //check if not already stored
 783         if(iprimary == tempo[kndex]){
 784           notIn = kFALSE ;
 785         }
 786       }
 787       if(notIn){
 788         if(fMulTrack<fMaxTrack){
 789           tempo[fMulTrack]=iprimary ;
 790           fMulTrack++ ;
 791         }
 792         else{
 793           Error("EvalPrimaries", "GetNprimaries ERROR > increase fMaxTrack!!!" ) ;
 794           break ;
 795         }
 796       }
 797     }
 798   } // all digits
 799   if(fMulTrack > 0){
 800     if(fTracksList)delete [] fTracksList;
 801     fTracksList = new Int_t[fMulTrack] ;
 802   }
 803   for(Int_t index = 0; index < fMulTrack; index++){
 804     fTracksList[index] = tempo[index] ;
 805   }
 806
 807   delete [] tempo ;
 808
 809 }
 810
 811 //____________________________________________________________________________
 812 void AliPHOSEmcRecPoint::EvalAll(TClonesArray * digits )
 813 {
 814 //   EvalCoreEnergy(logWeight, digits);
 815   EvalTime(digits) ;
 816   EvalPrimaries(digits) ;
 817   AliPHOSRecPoint::EvalAll(digits) ;
 818 }
 819 //____________________________________________________________________________
 820 void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TVector3 &vtx, TClonesArray * digits )
 821 {
 822   // Evaluates all shower parameters
 823   TVector3 vInc ;
 824   EvalLocalPosition(logWeight, vtx, digits,vInc) ;
 825   EvalElipsAxis(logWeight, digits, vInc) ; //they are evaluated with momenta
 826   EvalMoments(logWeight, digits, vInc) ;
 827   EvalDispersion(logWeight, digits, vInc) ;
 828 }
 829 //____________________________________________________________________________
 830 void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TVector3 &vtx, TClonesArray * digits, TVector3 &vInc)
 831 {
 832   // Calculates the center of gravity in the local PHOS-module coordinates
 833   Float_t wtot = 0. ;
 834
 835   Int_t relid[4] ;
 836
 837   Float_t x = 0. ;
 838   Float_t z = 0. ;
 839
 840   AliPHOSDigit * digit ;
 841
 842   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 843
 844   Int_t iDigit;
 845
 846   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 847     digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ;
 848
 849     Float_t xi ;
 850     Float_t zi ;
 851     phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 852     phosgeom->RelPosInModule(relid, xi, zi);
 853     if (fAmp>0 && fEnergyList[iDigit]>0) {
 854       Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
 855       x    += xi * w ;
 856       z    += zi * w ;
 857       wtot += w ;
 858     }
 859     else
 860       AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp));
 861   }
 862   if (wtot>0) {
 863     x /= wtot ;
 864     z /= wtot ;
 865   }
 866   else
 867     AliError(Form("Wrong weight %f\n", wtot));
 868
 869   // Correction for the depth of the shower starting point (TDR p 127)
 870   Float_t para = 0.925 ;
 871   Float_t parb = 6.52 ;
 872
 873   phosgeom->GetIncidentVector(vtx,GetPHOSMod(),x,z,vInc) ;
 874
 875   Float_t depthx = 0.;
 876   Float_t depthz = 0.;
 877   if (fAmp>0&&vInc.Y()!=0.) {
 878     depthx = ( para * TMath::Log(fAmp) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
 879     depthz = ( para * TMath::Log(fAmp) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
 880   }
 881   else
 882     AliError(Form("Wrong amplitude %f\n", fAmp));
 883
 884   fLocPos.SetX(x - depthx)  ;
 885   fLocPos.SetY(0.) ;
 886   fLocPos.SetZ(z - depthz)  ;
 887
 888 }
 889
 890 //____________________________________________________________________________
 891 Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) const
 892 {
 893   // Finds the maximum energy in the cluster
 894
 895   Float_t menergy = 0. ;
 896
 897   Int_t iDigit;
 898
 899   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 900
 901     if(fEnergyList[iDigit] > menergy)
 902       menergy = fEnergyList[iDigit] ;
 903   }
 904   return menergy ;
 905 }
 906
 907 //____________________________________________________________________________
 908 Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) const
 909 {
 910   // Calculates the multiplicity of digits with energy larger than H*energy
 911
 912   Int_t multipl   = 0 ;
 913   Int_t iDigit ;
 914   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 915
 916     if(fEnergyList[iDigit] > H * fAmp)
 917       multipl++ ;
 918   }
 919   return multipl ;
 920 }
 921
 922 //____________________________________________________________________________
 923 Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax( AliPHOSDigit **  maxAt, Float_t * maxAtEnergy,
 924                                                Float_t locMaxCut,TClonesArray * digits) const
 925 {
 926   // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
 927   // energy difference between two local maxima
 928
 929   AliPHOSDigit * digit ;
 930   AliPHOSDigit * digitN ;
 931
 932
 933   Int_t iDigitN ;
 934   Int_t iDigit ;
 935
 936   for(iDigit = 0; iDigit < fMulDigit; iDigit++)
 937     maxAt[iDigit] = (AliPHOSDigit*) digits->At(fDigitsList[iDigit])  ;
 938
 939
 940   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
 941     if(maxAt[iDigit]) {
 942       digit = maxAt[iDigit] ;
 943
 944       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
 945         if(iDigit == iDigitN)
 946           continue ;
 947
 948         digitN = (AliPHOSDigit *) digits->At(fDigitsList[iDigitN]) ;
 949
 950         if ( AreNeighbours(digit, digitN) ) {
 951           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
 952             maxAt[iDigitN] = 0 ;
 953             // but may be digit too is not local max ?
 954             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut)
 955               maxAt[iDigit] = 0 ;
 956           }
 957           else {
 958             maxAt[iDigit] = 0 ;
 959             // but may be digitN too is not local max ?
 960             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut)
 961               maxAt[iDigitN] = 0 ;
 962           }
 963         } // if Areneighbours
 964       } // while digitN
 965     } // slot not empty
 966   } // while digit
 967
 968   iDigitN = 0 ;
 969   for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
 970     if(maxAt[iDigit]){
 971       maxAt[iDigitN] = maxAt[iDigit] ;
 972       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
 973       iDigitN++ ;
 974     }
 975   }
 976   return iDigitN ;
 977 }
 978 //____________________________________________________________________________
 979 void AliPHOSEmcRecPoint::EvalTime(TClonesArray * /*digits*/)
 980 {
 981   // Define a rec.point time as a time in the cell with the maximum energy
 982   //Time already evaluated during AddDigit()
 983
 984 /*
 985   Float_t maxE = 0;
 986   Int_t maxAt = 0;
 987   for(Int_t idig=0; idig < fMulDigit; idig++){
 988     if(fEnergyList[idig] > maxE){
 989       maxE = fEnergyList[idig] ;
 990       maxAt = idig;
 991     }
 992   }
 993   fTime = ((AliPHOSDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
 994 */
 995 }
 996 //____________________________________________________________________________
 997 void AliPHOSEmcRecPoint::Purify(Float_t threshold, const TClonesArray * digits){
 998   //Removes digits below threshold
 999
1000   Int_t tempo[fMaxDigit];
1001   Float_t tempoE[fMaxDigit];
1002
1003   Int_t mult = 0 ;
1004   for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){
1005     if(fEnergyList[iDigit] > threshold){
1006       tempo[mult] = fDigitsList[iDigit] ;
1007       tempoE[mult] = fEnergyList[iDigit] ;
1008       mult++ ;
1009     }
1010   }
1011
1012   if(mult==0){ //too soft cluster
1013     fMulDigit =0 ;
1014     fAmp = 0.; //Recalculate total energy
1015   }
1016
1017   //Remove non-connected cells
1018   Int_t index[mult] ;
1019   Bool_t used[mult] ;
1020   for(Int_t i=0; i<mult; i++) used[i]=0 ;
1021   Int_t inClu=0 ;
1022   Double_t eMax=0. ;
1023   //find maximum
1024   for(Int_t iDigit=0; iDigit<mult; iDigit++) {
1025     if(eMax<tempoE[iDigit]){
1026       eMax=tempoE[iDigit];
1027       index[0]=iDigit ;
1028       inClu=1 ;
1029     }
1030   }
1031   if(mult>0)
1032     used[index[0]]=kTRUE ; //mark as used
1033   for(Int_t i=0; i<inClu; i++){
1034     AliPHOSDigit * digit = (AliPHOSDigit *) digits->At(tempo[index[i]]) ;
1035     for(Int_t iDigit=0 ;iDigit<mult; iDigit++){
1036       if(used[iDigit]) //already used
1037         continue ;
1038       AliPHOSDigit * digitN = (AliPHOSDigit *) digits->At(tempo[iDigit]) ;
1039       if(AreNeighbours(digit,digitN)){
1040         index[inClu]= iDigit ;
1041         inClu++ ;
1042         used[iDigit]=kTRUE ;
1043       }
1044     }
1045   }
1046
1047   fMulDigit = inClu ;
1048   delete [] fDigitsList ;
1049   delete [] fEnergyList ;
1050   fDigitsList = new Int_t[fMulDigit];
1051   fEnergyList = new Float_t[fMulDigit];
1052
1053   fAmp = 0.; //Recalculate total energy
1054   for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){
1055      fDigitsList[iDigit] = tempo[index[iDigit]];
1056      fEnergyList[iDigit] = tempoE[index[iDigit]] ;
1057      fAmp+=tempoE[index[iDigit]];
1058   }
1059 }
1060 //____________________________________________________________________________
1061 void AliPHOSEmcRecPoint::Print(Option_t *) const
1062 {
1063   // Print the list of digits belonging to the cluster
1064
1065   TString message ;
1066   message  = "AliPHOSEmcRecPoint:\n" ;
1067   message += "Digit multiplicity = %d" ;
1068   message += ", cluster Energy = %f" ;
1069   message += ", number of primaries = %d" ;
1070   message += ", rec.point index = %d \n" ;
1071   printf(message.Data(), fMulDigit, fAmp, fMulTrack,GetIndexInList() ) ;
1072
1073   Int_t iDigit;
1074   printf(" digits ids = ") ;
1075   for(iDigit=0; iDigit<fMulDigit; iDigit++)
1076     printf(" %d ", fDigitsList[iDigit] ) ;
1077
1078   printf("\n digit energies = ") ;
1079   for(iDigit=0; iDigit<fMulDigit; iDigit++)
1080     printf(" %f ", fEnergyList[iDigit] ) ;
1081
1082   printf("\n digit primaries  ") ;
1083   if (fMulTrack<1) printf("... no primaries");
1084   for(iDigit = 0;iDigit < fMulTrack; iDigit++)
1085     printf(" %d ", fTracksList[iDigit]) ;
1086   printf("\n") ;
1087
1088 }
1089
1090