EMCAL/AliEMCALRecPoint.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 /* $Id$ */
  16 //_________________________________________________________________________
  17 //  Reconstructed Points for the EMCAL
  18 //  A RecPoint is a cluster of digits
  19 //*-- Author: Yves Schutz (SUBATECH)
  20 //*-- Author: Dmitri Peressounko (RRC KI & SUBATECH)
  21 //*-- Author: Heather Gray (LBL) merged AliEMCALRecPoint and AliEMCALTowerRecPoint 02/04
  22
  23 // --- ROOT system ---
  24 class Riostream;
  25 #include <TPad.h>
  26 class TGraph;
  27 class TPaveText;
  28 #include <TClonesArray.h>
  29 #include <TMath.h>
  30
  31 // --- Standard library ---
  32
  33 // --- AliRoot header files ---
  34 //#include "AliGenerator.h"
  35 class AliGenerator;
  36 #include "AliRunLoader.h"
  37 #include "AliRun.h"
  38 class AliEMCAL;
  39 #include "AliEMCALLoader.h"
  40 #include "AliEMCALGeometry.h"
  41 #include "AliEMCALHit.h"
  42 #include "AliEMCALDigit.h"
  43 #include "AliEMCALRecPoint.h"
  44
  45 ClassImp(AliEMCALRecPoint)
  46
  47 //____________________________________________________________________________
  48 AliEMCALRecPoint::AliEMCALRecPoint()
  49   : AliRecPoint()
  50 {
  51   // ctor
  52   fClusterType = -1;
  53   fMaxTrack = 0 ;
  54   fMulDigit   = 0 ;
  55   fMaxParent = 0;
  56   fMulParent = 0;
  57   fAmp   = 0. ;
  58   fCoreEnergy = 0 ;
  59   fEnergyList = 0 ;
  60   fTimeList = 0 ;
  61   fAbsIdList  = 0;
  62   fParentsList = 0;
  63   fTime = 0. ;
  64   //  fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
  65   fCoreRadius = 10;        //HG Check this
  66
  67   AliRunLoader *rl = AliRunLoader::GetRunLoader();
  68   fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
  69   //fGeomPtr = AliEMCALGeometry::GetInstance();
  70   fGeomPtr->GetTransformationForSM(); // Global <-> Local
  71 }
  72
  73 //____________________________________________________________________________
  74 AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) : AliRecPoint(opt)
  75 {
  76   // ctor
  77   fClusterType = -1;
  78   fMaxTrack = 1000 ;
  79   fMaxParent = 1000;
  80   fMulDigit   = 0 ;
  81   fMulParent = 0;
  82   fAmp   = 0. ;
  83   fCoreEnergy = 0 ;
  84   fEnergyList = 0 ;
  85   fTimeList = 0 ;
  86   fAbsIdList  = 0;
  87   fParentsList = new Int_t[fMaxParent];
  88   fTime = -1. ;
  89   //fLocPos.SetX(1.e+6)  ;      //Local position should be evaluated
  90   fCoreRadius = 10;        //HG Check this
  91   //fGeomPtr = AliEMCALGeometry::GetInstance();
  92   AliRunLoader *rl = AliRunLoader::GetRunLoader();
  93   fGeomPtr = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
  94   fGeomPtr->GetTransformationForSM(); // Global <-> Local
  95 }
  96 //____________________________________________________________________________
  97 AliEMCALRecPoint::~AliEMCALRecPoint()
  98 {
  99   // dtor
 100   if ( fEnergyList )
 101     delete[] fEnergyList ;
 102   if ( fTimeList )
 103     delete[] fTimeList ;
 104   if ( fAbsIdList )
 105     delete[] fAbsIdList ;
 106    if ( fParentsList)
 107     delete[] fParentsList;
 108 }
 109
 110 //____________________________________________________________________________
 111 void AliEMCALRecPoint::AddDigit(AliEMCALDigit & digit, Float_t Energy)
 112 {
 113   // Adds a digit to the RecPoint
 114   // and accumulates the total amplitude and the multiplicity
 115
 116   if(fEnergyList == 0)
 117     fEnergyList =  new Float_t[fMaxDigit];
 118   if(fTimeList == 0)
 119     fTimeList =  new Float_t[fMaxDigit];
 120   if(fAbsIdList == 0) {
 121     fAbsIdList =  new Int_t[fMaxDigit];
 122     fSuperModuleNumber = fGeomPtr->GetSuperModuleNumber(digit.GetId());
 123   }
 124
 125   if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
 126     fMaxDigit*=2 ;
 127     Int_t * tempo = new Int_t[fMaxDigit];
 128     Float_t * tempoE =  new Float_t[fMaxDigit];
 129     Float_t * tempoT =  new Float_t[fMaxDigit];
 130     Int_t * tempoId = new Int_t[fMaxDigit];
 131
 132     Int_t index ;
 133     for ( index = 0 ; index < fMulDigit ; index++ ){
 134       tempo[index]   = fDigitsList[index] ;
 135       tempoE[index]  = fEnergyList[index] ;
 136       tempoT[index]  = fTimeList[index] ;
 137       tempoId[index] = fAbsIdList[index] ;
 138     }
 139
 140     delete [] fDigitsList ;
 141     fDigitsList =  new Int_t[fMaxDigit];
 142
 143     delete [] fEnergyList ;
 144     fEnergyList =  new Float_t[fMaxDigit];
 145
 146     delete [] fTimeList ;
 147     fTimeList =  new Float_t[fMaxDigit];
 148
 149     delete [] fAbsIdList ;
 150     fAbsIdList =  new Int_t[fMaxDigit];
 151
 152     for ( index = 0 ; index < fMulDigit ; index++ ){
 153       fDigitsList[index] = tempo[index] ;
 154       fEnergyList[index] = tempoE[index] ;
 155       fTimeList[index] = tempoT[index] ;
 156       fAbsIdList[index]  = tempoId[index] ;
 157     }
 158
 159     delete [] tempo ;
 160     delete [] tempoE ;
 161     delete [] tempoT ;
 162     delete [] tempoId ;
 163   } // if
 164
 165   fDigitsList[fMulDigit]   = digit.GetIndexInList()  ;
 166   fEnergyList[fMulDigit]   = Energy ;
 167   fTimeList[fMulDigit]     = digit.GetTime() ;
 168   fAbsIdList[fMulDigit]    = digit.GetId();
 169   fMulDigit++ ;
 170   fAmp += Energy ;
 171
 172 }
 173 //____________________________________________________________________________
 174 Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * digit2 ) const
 175 {
 176   // Tells if (true) or not (false) two digits are neighbours
 177   // A neighbour is defined as being two digits which share a corner
 178
 179   static Bool_t areNeighbours = kFALSE ;
 180   static Int_t nSupMod=0, nTower=0, nIphi=0, nIeta=0;
 181   static int nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0;
 182   static Int_t relid1[2] , relid2[2] ; // ieta, iphi
 183   static Int_t rowdiff=0, coldiff=0;
 184
 185   areNeighbours = kFALSE ;
 186
 187   fGeomPtr->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta);
 188   fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, relid1[0],relid1[1]);
 189
 190   fGeomPtr->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1);
 191   fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, relid2[0],relid2[1]);
 192
 193   rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ;
 194   coldiff = TMath::Abs( relid1[1] - relid2[1] ) ;
 195
 196   if (( coldiff <= 1 )  && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
 197   areNeighbours = kTRUE ;
 198
 199   return areNeighbours;
 200 }
 201
 202 //____________________________________________________________________________
 203 Int_t AliEMCALRecPoint::Compare(const TObject * obj) const
 204 {
 205   // Compares two RecPoints according to their position in the EMCAL modules
 206
 207   Float_t delta = 1 ; //Width of "Sorting row". If you change this
 208                       //value (what is senseless) change as well delta in
 209                       //AliEMCALTrackSegmentMakerv* and other RecPoints...
 210   Int_t rv ;
 211
 212   AliEMCALRecPoint * clu = (AliEMCALRecPoint *)obj ;
 213
 214   TVector3 locpos1;
 215   GetLocalPosition(locpos1);
 216   TVector3 locpos2;
 217   clu->GetLocalPosition(locpos2);
 218
 219   Int_t rowdif = (Int_t)(TMath::Ceil(locpos1.X()/delta)-TMath::Ceil(locpos2.X()/delta)) ;
 220   if (rowdif> 0)
 221     rv = 1 ;
 222   else if(rowdif < 0)
 223     rv = -1 ;
 224   else if(locpos1.Y()>locpos2.Y())
 225     rv = -1 ;
 226   else
 227     rv = 1 ;
 228
 229   return rv ;
 230 }
 231
 232 //____________________________________________________________________________
 233 Int_t AliEMCALRecPoint::DistancetoPrimitive(Int_t px, Int_t py)
 234 {
 235   // Compute distance from point px,py to  a AliEMCALRecPoint considered as a Tmarker
 236   // Compute the closest distance of approach from point px,py to this marker.
 237   // The distance is computed in pixels units.
 238   // HG Still need to update -> Not sure what this should achieve
 239
 240   TVector3 pos(0.,0.,0.) ;
 241   GetLocalPosition(pos) ;
 242   Float_t x =  pos.X() ;
 243   Float_t y =  pos.Y() ;
 244   const Int_t kMaxDiff = 10;
 245   Int_t pxm  = gPad->XtoAbsPixel(x);
 246   Int_t pym  = gPad->YtoAbsPixel(y);
 247   Int_t dist = (px-pxm)*(px-pxm) + (py-pym)*(py-pym);
 248
 249   if (dist > kMaxDiff) return 9999;
 250   return dist;
 251 }
 252
 253 //___________________________________________________________________________
 254  void AliEMCALRecPoint::Draw(Option_t *option)
 255  {
 256    // Draw this AliEMCALRecPoint with its current attributes
 257
 258    AppendPad(option);
 259  }
 260
 261 //______________________________________________________________________________
 262 void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t)
 263 {
 264   // Execute action corresponding to one event
 265   // This member function is called when a AliEMCALRecPoint is clicked with the locator
 266   //
 267   // If Left button is clicked on AliEMCALRecPoint, the digits are switched on
 268   // and switched off when the mouse button is released.
 269
 270   //  static Int_t pxold, pyold;
 271
 272   /* static TGraph *  digitgraph = 0 ;
 273   static TPaveText* clustertext = 0 ;
 274
 275   if (!gPad->IsEditable()) return;
 276
 277   switch (event) {
 278
 279
 280   case kButton1Down:{
 281     AliEMCALDigit * digit ;
 282     AliEMCALGeometry * emcalgeom =  (AliEMCALGetter::Instance())->EMCALGeometry() ;
 283
 284     Int_t iDigit;
 285     Int_t relid[2] ;
 286
 287     const Int_t kMulDigit=AliEMCALRecPoint::GetDigitsMultiplicity() ;
 288     Float_t * xi = new Float_t [kMulDigit] ;
 289     Float_t * zi = new Float_t [kMulDigit] ;
 290
 291     for(iDigit = 0; iDigit < kMulDigit; iDigit++) {
 292       Fatal("AliEMCALRecPoint::ExecuteEvent", " -> Something wrong with the code");
 293       digit = 0 ; //dynamic_cast<AliEMCALDigit *>((fDigitsList)[iDigit]);
 294       emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ;
 295       emcalgeom->PosInAlice(relid, xi[iDigit], zi[iDigit]) ;
 296     }
 297
 298     if (!digitgraph) {
 299       digitgraph = new TGraph(fMulDigit,xi,zi);
 300       digitgraph-> SetMarkerStyle(5) ;
 301       digitgraph-> SetMarkerSize(1.) ;
 302       digitgraph-> SetMarkerColor(1) ;
 303       digitgraph-> Draw("P") ;
 304     }
 305     if (!clustertext) {
 306
 307       TVector3 pos(0.,0.,0.) ;
 308       GetLocalPosition(pos) ;
 309       clustertext = new TPaveText(pos.X()-10,pos.Z()+10,pos.X()+50,pos.Z()+35,"") ;
 310       Text_t  line1[40] ;
 311       Text_t  line2[40] ;
 312       sprintf(line1,"Energy=%1.2f GeV",GetEnergy()) ;
 313       sprintf(line2,"%d Digits",GetDigitsMultiplicity()) ;
 314       clustertext ->AddText(line1) ;
 315       clustertext ->AddText(line2) ;
 316       clustertext ->Draw("");
 317     }
 318     gPad->Update() ;
 319     Print("") ;
 320     delete[] xi ;
 321     delete[] zi ;
 322    }
 323
 324 break;
 325
 326   case kButton1Up:
 327     if (digitgraph) {
 328       delete digitgraph  ;
 329       digitgraph = 0 ;
 330     }
 331     if (clustertext) {
 332       delete clustertext ;
 333       clustertext = 0 ;
 334     }
 335
 336     break;
 337
 338     }*/
 339 }
 340 //____________________________________________________________________________
 341 void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits)
 342 {
 343   // Evaluates all shower parameters
 344
 345   EvalLocalPosition(logWeight, digits) ;
 346   //  printf("eval position done\n");
 347   EvalElipsAxis(logWeight, digits) ;
 348   //  printf("eval axis done\n");
 349   EvalDispersion(logWeight, digits) ;
 350   //  printf("eval dispersion done\n");
 351   //EvalCoreEnergy(logWeight, digits);
 352  // printf("eval energy done\n");
 353   EvalTime(digits) ;
 354   //  printf("eval time done\n");
 355
 356   EvalPrimaries(digits) ;
 357   //  printf("eval pri done\n");
 358   EvalParents(digits);
 359   //  printf("eval parent done\n");
 360 }
 361
 362 //____________________________________________________________________________
 363 void  AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits)
 364 {
 365   // Calculates the dispersion of the shower at the origin of the RecPoint
 366
 367   Double_t d = 0., wtot = 0., w = 0., xyzi[3], diff=0.;
 368   Int_t iDigit=0, nstat=0, i=0;
 369   AliEMCALDigit * digit ;
 370
 371   // Calculates the centre of gravity in the local EMCAL-module coordinates
 372   if (!fLocPos.Mag())
 373     EvalLocalPosition(logWeight, digits) ;
 374
 375   // Calculates the dispersion in coordinates
 376   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
 377     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
 378
 379     fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
 380     w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 381
 382     if(w>0.0) {
 383       wtot += w;
 384       nstat++;
 385       for(i=0; i<3; i++ ) {
 386         diff = xyzi[i] - double(fLocPos[i]);
 387         d   += w * diff*diff;
 388       }
 389     }
 390   }
 391
 392   if ( wtot > 0 && nstat>1) d /= wtot ;
 393   else                      d = 0. ;
 394
 395   fDispersion = TMath::Sqrt(d) ;
 396 }
 397
 398 //____________________________________________________________________________
 399 void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits)
 400 {
 401   // Calculates the center of gravity in the local EMCAL-module coordinates
 402   //  Info("Print", " logWeight %f : cluster energy %f ", logWeight, fAmp); // for testing
 403
 404   AliEMCALDigit * digit;
 405   Int_t i=0, nstat=0;
 406   Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
 407
 408   for(Int_t iDigit=0; iDigit<fMulDigit; iDigit++) {
 409     digit = dynamic_cast<AliEMCALDigit *>(digits->At(fDigitsList[iDigit])) ;
 410
 411     fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
 412     // printf(" Id %i : Local x,y,z %f %f %f \n", digit->GetId(), xyzi[0], xyzi[1], xyzi[2]);
 413
 414     if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ));
 415     else                w = fEnergyList[iDigit]; // just energy
 416
 417     if(w>0.0) {
 418       wtot += w ;
 419       nstat++;
 420       for(i=0; i<3; i++ ) {
 421         clXYZ[i]    += (w*xyzi[i]);
 422         clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]);
 423       }
 424     }
 425   }
 426   //  cout << " wtot " << wtot << endl;
 427   if ( wtot > 0 ) {
 428     //    xRMS   = TMath::Sqrt(x2m - xMean*xMean);
 429     for(i=0; i<3; i++ ) {
 430       clXYZ[i] /= wtot;
 431       if(nstat>1) {
 432         clRmsXYZ[i] /= (wtot*wtot);
 433         clRmsXYZ[i]  =  clRmsXYZ[i] - clXYZ[i]*clXYZ[i];
 434         if(clRmsXYZ[i] > 0.0) {
 435           clRmsXYZ[i] =  TMath::Sqrt(clRmsXYZ[i]);
 436         } else      clRmsXYZ[i] = 0;
 437       } else        clRmsXYZ[i] = 0;
 438     }
 439   } else {
 440     for(i=0; i<3; i++ ) {
 441       clXYZ[i] = clRmsXYZ[i] = -1.;
 442     }
 443   }
 444   // clRmsXYZ[i] ??
 445   fLocPos.SetX(clXYZ[0]);
 446   fLocPos.SetY(clXYZ[1]);
 447   fLocPos.SetZ(clXYZ[2]);
 448
 449 //  if (gDebug==2)
 450 //    printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ;
 451   fLocPosM = 0 ; // covariance matrix
 452 }
 453
 454 //void AliEMCALRecPoint::EvalLocalPositionSimple()
 455 //{ // Weight is proportional of cell energy
 456 //}
 457
 458 //______________________________________________________________________________
 459 void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits)
 460 {
 461   // This function calculates energy in the core,
 462   // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius
 463   // in accordance with shower profile the energy deposition
 464   // should be less than 2%
 465
 466   AliEMCALDigit * digit ;
 467   const Float_t kDeg2Rad = TMath::DegToRad() ;
 468
 469   Int_t iDigit;
 470
 471   if (!fLocPos.Mag()) {
 472     EvalLocalPosition(logWeight, digits);
 473   }
 474
 475   for(iDigit=0; iDigit < fMulDigit; iDigit++) {
 476     digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ;
 477
 478     Float_t ieta = 0.0;
 479     Float_t iphi = 0.0;
 480     //fGeomPtr->PosInAlice(digit->GetId(), ieta, iphi);
 481     fGeomPtr->EtaPhiFromIndex(digit->GetId(),ieta, iphi) ;
 482     iphi = iphi * kDeg2Rad;
 483
 484     Float_t distance = TMath::Sqrt((ieta-fLocPos.X())*(ieta-fLocPos.X())+(iphi-fLocPos.Y())*(iphi-fLocPos.Y())) ;
 485     if(distance < fCoreRadius)
 486       fCoreEnergy += fEnergyList[iDigit] ;
 487   }
 488
 489 }
 490 //____________________________________________________________________________
 491 void  AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits)
 492 {
 493   // Calculates the axis of the shower ellipsoid in eta and phi
 494
 495   Double_t wtot = 0. ;
 496   Double_t x    = 0.;
 497   Double_t z    = 0.;
 498   Double_t dxx  = 0.;
 499   Double_t dzz  = 0.;
 500   Double_t dxz  = 0.;
 501
 502   const Float_t kDeg2Rad = TMath::DegToRad();
 503   AliEMCALDigit * digit ;
 504
 505   TString gn(fGeomPtr->GetName());
 506
 507   Int_t iDigit;
 508
 509   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 510     digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit])  ;
 511     Float_t etai = 0. ;
 512     Float_t phii = 0. ;
 513     if(gn.Contains("SHISH")) { // have to be change - Feb 28, 2006
 514     //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion
 515     // for this geometry does not exist
 516       int nSupMod=0, nTower=0, nIphi=0, nIeta=0;
 517       int iphi=0, ieta=0;
 518       fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta);
 519       fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta);
 520       etai=(Float_t)ieta;
 521       phii=(Float_t)iphi;
 522     } else {
 523       fGeomPtr->EtaPhiFromIndex(digit->GetId(), etai, phii);
 524       phii = phii * kDeg2Rad;
 525     }
 526
 527     Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
 528
 529     dxx  += w * etai * etai ;
 530     x    += w * etai ;
 531     dzz  += w * phii * phii ;
 532     z    += w * phii ;
 533
 534     dxz  += w * etai * phii ;
 535
 536     wtot += w ;
 537   }
 538
 539   if ( wtot > 0 ) {
 540     dxx /= wtot ;
 541     x   /= wtot ;
 542     dxx -= x * x ;
 543     dzz /= wtot ;
 544     z   /= wtot ;
 545     dzz -= z * z ;
 546     dxz /= wtot ;
 547     dxz -= x * z ;
 548
 549     fLambda[0] =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 550     if(fLambda[0] > 0)
 551       fLambda[0] = TMath::Sqrt(fLambda[0]) ;
 552     else
 553       fLambda[0] = 0;
 554
 555     fLambda[1] =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 556
 557     if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda.
 558       fLambda[1] = TMath::Sqrt(fLambda[1]) ;
 559     else
 560       fLambda[1]= 0. ;
 561   } else {
 562     fLambda[0]= 0. ;
 563     fLambda[1]= 0. ;
 564   }
 565
 566   //  printf("Evalaxis: lambdas  = %f,%f", fLambda[0],fLambda[1]) ;
 567
 568 }
 569
 570 //______________________________________________________________________________
 571 void  AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits)
 572 {
 573   // Constructs the list of primary particles (tracks) which have contributed to this RecPoint
 574
 575   AliEMCALDigit * digit ;
 576   Int_t * tempo    = new Int_t[fMaxTrack] ;
 577
 578   Int_t index ;
 579   for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
 580     digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ;
 581     Int_t nprimaries = digit->GetNprimary() ;
 582     if ( nprimaries == 0 ) continue ;
 583     Int_t * newprimaryarray = new Int_t[nprimaries] ;
 584     Int_t ii ;
 585     for ( ii = 0 ; ii < nprimaries ; ii++)
 586       newprimaryarray[ii] = digit->GetPrimary(ii+1) ;
 587
 588     Int_t jndex ;
 589     for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit
 590       if ( fMulTrack > fMaxTrack ) {
 591         fMulTrack = fMaxTrack ;
 592         Error("GetNprimaries", "increase fMaxTrack ")  ;
 593         break ;
 594       }
 595       Int_t newprimary = newprimaryarray[jndex] ;
 596       Int_t kndex ;
 597       Bool_t already = kFALSE ;
 598       for ( kndex = 0 ; kndex < fMulTrack ; kndex++ ) { //check if not already stored
 599         if ( newprimary == tempo[kndex] ){
 600           already = kTRUE ;
 601           break ;
 602         }
 603       } // end of check
 604       if ( !already && (fMulTrack < fMaxTrack)) { // store it
 605         tempo[fMulTrack] = newprimary ;
 606         fMulTrack++ ;
 607       } // store it
 608     } // all primaries in digit
 609     delete [] newprimaryarray ;
 610   } // all digits
 611
 612
 613   fTracksList = new Int_t[fMulTrack] ;
 614   for(index = 0; index < fMulTrack; index++)
 615    fTracksList[index] = tempo[index] ;
 616
 617   delete [] tempo ;
 618
 619 }
 620
 621 //______________________________________________________________________________
 622 void  AliEMCALRecPoint::EvalParents(TClonesArray * digits)
 623 {
 624   // Constructs the list of parent particles (tracks) which have contributed to this RecPoint
 625
 626   AliEMCALDigit * digit ;
 627   Int_t * tempo    = new Int_t[fMaxParent] ;
 628
 629   Int_t index ;
 630   for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits
 631     digit = dynamic_cast<AliEMCALDigit *>(digits->At( fDigitsList[index] )) ;
 632     Int_t nparents = digit->GetNiparent() ;
 633     if ( nparents == 0 ) continue ;
 634     Int_t * newparentarray = new Int_t[nparents] ;
 635     Int_t ii ;
 636     for ( ii = 0 ; ii < nparents ; ii++)
 637       newparentarray[ii] = digit->GetIparent(ii+1) ;
 638
 639     Int_t jndex ;
 640     for ( jndex = 0 ; jndex < nparents ; jndex++ ) { // all primaries in digit
 641       if ( fMulParent > fMaxParent ) {
 642         fMulTrack = - 1 ;
 643         Error("GetNiparent", "increase fMaxParent")  ;
 644         break ;
 645       }
 646       Int_t newparent = newparentarray[jndex] ;
 647       Int_t kndex ;
 648       Bool_t already = kFALSE ;
 649       for ( kndex = 0 ; kndex < fMulParent ; kndex++ ) { //check if not already stored
 650         if ( newparent == tempo[kndex] ){
 651           already = kTRUE ;
 652           break ;
 653         }
 654       } // end of check
 655       if ( !already && (fMulTrack < fMaxTrack)) { // store it
 656         tempo[fMulParent] = newparent ;
 657         fMulParent++ ;
 658       } // store it
 659     } // all parents in digit
 660     delete [] newparentarray ;
 661   } // all digits
 662
 663   if (fMulParent>0) {
 664     fParentsList = new Int_t[fMulParent] ;
 665     for(index = 0; index < fMulParent; index++)
 666       fParentsList[index] = tempo[index] ;
 667   }
 668
 669   delete [] tempo ;
 670
 671 }
 672
 673 //____________________________________________________________________________
 674 void AliEMCALRecPoint::GetLocalPosition(TVector3 & lpos) const
 675 {
 676   // returns the position of the cluster in the local reference system of ALICE
 677
 678   lpos.SetX(fLocPos.X()) ;
 679   lpos.SetY(fLocPos.Y()) ;
 680   lpos.SetZ(fLocPos.Z()) ;
 681 }
 682
 683 //____________________________________________________________________________
 684 void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const
 685 {
 686   // returns the position of the cluster in the global reference system of ALICE
 687   // These are now the Cartesian X, Y and Z
 688   //  cout<<" geom "<<geom<<endl;
 689   fGeomPtr->GetGlobal(fLocPos, gpos, fSuperModuleNumber);
 690 }
 691
 692 //____________________________________________________________________________
 693 Float_t AliEMCALRecPoint::GetMaximalEnergy(void) const
 694 {
 695   // Finds the maximum energy in the cluster
 696
 697   Float_t menergy = 0. ;
 698
 699   Int_t iDigit;
 700
 701   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 702
 703     if(fEnergyList[iDigit] > menergy)
 704       menergy = fEnergyList[iDigit] ;
 705   }
 706   return menergy ;
 707 }
 708
 709 //____________________________________________________________________________
 710 Int_t AliEMCALRecPoint::GetMultiplicityAtLevel(Float_t H) const
 711 {
 712   // Calculates the multiplicity of digits with energy larger than H*energy
 713
 714   Int_t multipl   = 0 ;
 715   Int_t iDigit ;
 716   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 717
 718     if(fEnergyList[iDigit] > H * fAmp)
 719       multipl++ ;
 720   }
 721   return multipl ;
 722 }
 723
 724 //____________________________________________________________________________
 725 Int_t  AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit **  maxAt, Float_t * maxAtEnergy,
 726                                                Float_t locMaxCut,TClonesArray * digits) const
 727 {
 728   // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum
 729   // energy difference between two local maxima
 730
 731   AliEMCALDigit * digit ;
 732   AliEMCALDigit * digitN ;
 733
 734   Int_t iDigitN ;
 735   Int_t iDigit ;
 736
 737   for(iDigit = 0; iDigit < fMulDigit; iDigit++)
 738     maxAt[iDigit] = (AliEMCALDigit*) digits->At(fDigitsList[iDigit])  ;
 739
 740   for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
 741     if(maxAt[iDigit]) {
 742       digit = maxAt[iDigit] ;
 743
 744       for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
 745         digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ;
 746
 747         if ( AreNeighbours(digit, digitN) ) {
 748           if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
 749             maxAt[iDigitN] = 0 ;
 750             // but may be digit too is not local max ?
 751             if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut)
 752               maxAt[iDigit] = 0 ;
 753           }
 754           else {
 755             maxAt[iDigit] = 0 ;
 756             // but may be digitN too is not local max ?
 757             if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut)
 758               maxAt[iDigitN] = 0 ;
 759           }
 760         } // if Areneighbours
 761       } // while digitN
 762     } // slot not empty
 763   } // while digit
 764
 765   iDigitN = 0 ;
 766   for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
 767     if(maxAt[iDigit] ){
 768       maxAt[iDigitN] = maxAt[iDigit] ;
 769       maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
 770       iDigitN++ ;
 771     }
 772   }
 773   return iDigitN ;
 774 }
 775
 776 //____________________________________________________________________________
 777 Int_t AliEMCALRecPoint::GetPrimaryIndex() const
 778 {
 779   // Get the primary track index in TreeK which deposits the most energy
 780   // in Digits which forms RecPoint. Kinematics, Hits and Digits must be
 781   // loaded before the call of the method.
 782
 783   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 784   if (!rl)
 785     AliError(Form(" No Runloader ")) ;
 786
 787   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
 788     (rl->GetDetectorLoader("EMCAL"));
 789
 790   // Get the list of digits forming this RecPoint
 791   Int_t  nDigits   = fMulDigit   ;
 792   Int_t *digitList = fDigitsList ;
 793
 794   // Find the digit with maximum amplitude
 795   AliEMCALDigit *digit = 0;
 796   TClonesArray *digits = emcalLoader->Digits();
 797   Int_t maxAmp = 0;
 798   Int_t bestDigitIndex = -1;
 799   for (Int_t iDigit=0; iDigit<nDigits; iDigit++) {
 800     digit = static_cast<AliEMCALDigit *>(digits->At(digitList[iDigit]));
 801     if (digit->GetAmp() > maxAmp) {
 802       maxAmp = digit->GetAmp();
 803       bestDigitIndex = iDigit;
 804     }
 805   }
 806
 807   digit = static_cast<AliEMCALDigit *>(digits->At(digitList[bestDigitIndex]));
 808
 809   // Get the list of hits producing this digit,
 810   // find which hit has deposited more energy
 811   // and find the primary track.
 812
 813   AliEMCALHit *hit = 0;
 814   TClonesArray *hits = emcalLoader->Hits();
 815
 816   Double_t maxedep  =  0;
 817   Int_t    maxtrack = -1;
 818   Int_t    nHits    = hits ->GetEntries();
 819   Int_t    id       = digit->GetId();
 820   for (Int_t iHit=0; iHit<nHits; iHit++) {
 821     hit = static_cast<AliEMCALHit*> (hits->At(iHit)) ;
 822     if(hit->GetId() == id){
 823       Double_t edep  = hit->GetEnergy();
 824       Int_t    track = hit->GetIparent();//Primary();
 825       if(edep > maxedep){
 826         maxedep  = edep;
 827         maxtrack = track;
 828       }
 829     }
 830   }
 831   if (maxtrack != -1) return maxtrack;
 832   return -12345;                       // no track found :(
 833 }
 834
 835 //____________________________________________________________________________
 836 void AliEMCALRecPoint::EvalTime(TClonesArray * digits){
 837   // time is set to the time of the digit with the maximum energy
 838
 839   Float_t maxE = 0;
 840   Int_t maxAt = 0;
 841   for(Int_t idig=0; idig < fMulDigit; idig++){
 842     if(fEnergyList[idig] > maxE){
 843       maxE = fEnergyList[idig] ;
 844       maxAt = idig;
 845     }
 846   }
 847   fTime = ((AliEMCALDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ;
 848
 849 }
 850
 851 //______________________________________________________________________________
 852 void AliEMCALRecPoint::Paint(Option_t *)
 853 {
 854   // Paint this ALiRecPoint as a TMarker  with its current attributes
 855
 856   TVector3 pos(0.,0.,0.)  ;
 857   GetLocalPosition(pos)   ;
 858   Coord_t x = pos.X()     ;
 859   Coord_t y = pos.Z()     ;
 860   Color_t markercolor = 1 ;
 861   Size_t  markersize = 1. ;
 862   Style_t markerstyle = 5 ;
 863
 864   if (!gPad->IsBatch()) {
 865     gVirtualX->SetMarkerColor(markercolor) ;
 866     gVirtualX->SetMarkerSize (markersize)  ;
 867     gVirtualX->SetMarkerStyle(markerstyle) ;
 868   }
 869   gPad->SetAttMarkerPS(markercolor,markerstyle,markersize) ;
 870   gPad->PaintPolyMarker(1,&x,&y,"") ;
 871 }
 872
 873 //______________________________________________________________________________
 874 Float_t AliEMCALRecPoint::EtaToTheta(Float_t arg) const
 875 {
 876   //Converts Theta (Radians) to Eta(Radians)
 877   return (2.*TMath::ATan(TMath::Exp(-arg)));
 878 }
 879
 880 //______________________________________________________________________________
 881 Float_t AliEMCALRecPoint::ThetaToEta(Float_t arg) const
 882 {
 883   //Converts Eta (Radians) to Theta(Radians)
 884   return (-1 * TMath::Log(TMath::Tan(0.5 * arg)));
 885 }
 886
 887 //____________________________________________________________________________
 888 void AliEMCALRecPoint::Print(Option_t *) const
 889 {
 890   // Print the list of digits belonging to the cluster
 891   return;
 892   TString message ;
 893   message  = "AliEMCALRecPoint:\n" ;
 894   message +=  " digits # = " ;
 895   Info("Print", message.Data()) ;
 896
 897   Int_t iDigit;
 898   for(iDigit=0; iDigit<fMulDigit; iDigit++)
 899     printf(" %d ", fDigitsList[iDigit] ) ;
 900   printf("\n");
 901
 902   Info("Print", " Energies = ") ;
 903   for(iDigit=0; iDigit<fMulDigit; iDigit++)
 904     printf(" %f ", fEnergyList[iDigit] ) ;
 905   printf("\n");
 906
 907   Info("Print", "\n Abs Ids = ") ;
 908   for(iDigit=0; iDigit<fMulDigit; iDigit++)
 909     printf(" %i ", fAbsIdList[iDigit] ) ;
 910   printf("\n");
 911
 912   Info("Print", " Primaries  ") ;
 913   for(iDigit = 0;iDigit < fMulTrack; iDigit++)
 914     printf(" %d ", fTracksList[iDigit]) ;
 915
 916   printf("\n Local x %6.2f y %7.2f z %7.1f \n", fLocPos[0], fLocPos[1], fLocPos[2]);
 917
 918   message  = "       ClusterType     = %d" ;
 919   message += "       Multiplicity    = %d" ;
 920   message += "       Cluster Energy  = %f" ;
 921   message += "       Core energy     = %f" ;
 922   message += "       Core radius     = %f" ;
 923   message += "       Number of primaries %d" ;
 924   message += "       Stored at position %d" ;
 925   Info("Print", message.Data(), fClusterType, fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ;
 926 }