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