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