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