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