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