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