[u/mrichter/AliRoot.git] / PHOS / AliPHOSEmcRecPoint.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//_________________________________________________________________________
// Rec Point in the PHOS EM calorimeter 
//*-- Author : Dmitri Peressounko RRC KI
//////////////////////////////////////////////////////////////////////////////

// --- ROOT system ---
#include "TPad.h"
#include "TH2.h"
#include "TMath.h" 
#include "TCanvas.h" 

// --- Standard library ---

#include <iostream>

// --- AliRoot header files ---

#include "AliPHOSGeometry.h"
#include "AliPHOSEmcRecPoint.h"
#include "AliRun.h"

ClassImp(AliPHOSEmcRecPoint)

//____________________________________________________________________________
AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
  : AliPHOSRecPoint()
{
  // ctor

  fMulDigit   = 0 ;  
  fAmp   = 0. ;   
  fEnergyList = new Float_t[fMaxDigit]; 
  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
  fDelta     =  PHOSGeom->GetCrystalSize(0) ; 
  fW0        = W0 ;          
  fLocMaxCut = LocMaxCut ; 
  fLocPos.SetX(1000000.)  ;      //Local position should be evaluated
}

//____________________________________________________________________________
AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() 
{
  // dtor 
}

//____________________________________________________________________________
void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
{
  // adds a digit to the digits list
  // and accumulates the total amplitude and the multiplicity 
  
  if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists 
    fMaxDigit*=2 ; 
    int * tempo = new ( int[fMaxDigit] ) ; 
    Float_t * tempoE =  new ( Float_t[fMaxDigit] ) ;

    Int_t index ;     
    for ( index = 0 ; index < fMulDigit ; index++ ){
      tempo[index] = fDigitsList[index] ;
      tempoE[index] = fEnergyList[index] ; 
    }
    
    delete [] fDigitsList ; 
    fDigitsList =  new ( int[fMaxDigit] ) ;
 
    delete [] fEnergyList ;
    fEnergyList =  new ( Float_t[fMaxDigit] ) ;

    for ( index = 0 ; index < fMulDigit ; index++ ){
      fDigitsList[index] = tempo[index] ;
      fEnergyList[index] = tempoE[index] ; 
    }
 
    delete [] tempo ;
    delete [] tempoE ; 
  } // if
  
  fDigitsList[fMulDigit]   =  (int) &digit  ; 
  fEnergyList[fMulDigit++] = Energy ;
  fAmp += Energy ; 
}

//____________________________________________________________________________
Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) 
{
  
  Bool_t aren = kFALSE ;
  
  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
  Int_t relid1[4] ; 
  PHOSGeom->AbsToRelNumbering(digit1->GetId(), relid1) ; 

  Int_t relid2[4] ; 
  PHOSGeom->AbsToRelNumbering(digit2->GetId(), relid2) ; 
  
  Int_t RowDiff = TMath::Abs( relid1[2] - relid2[2] ) ;  
  Int_t ColDiff = TMath::Abs( relid1[3] - relid2[3] ) ;  

  if (( ColDiff<=1 )  && ( RowDiff <= 1 ) && (ColDiff+RowDiff > 0)) 
    aren = kTRUE ;
  
  return aren ;
}

//____________________________________________________________________________
Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
{
  Int_t rv ; 

  AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; 

 
  Int_t PHOSMod1 = this->GetPHOSMod() ;
  Int_t PHOSMod2 = clu->GetPHOSMod() ;

  TVector3 LocPos1; 
  this->GetLocalPosition(LocPos1) ;
  TVector3 LocPos2;  
  clu->GetLocalPosition(LocPos2) ;  

  if(PHOSMod1 == PHOSMod2 ) {
    Int_t rowdif = (Int_t)TMath::Ceil(LocPos1.X()/fDelta)-(Int_t)TMath::Ceil(LocPos2.X()/fDelta) ;
    if (rowdif> 0) 
      rv = -1 ;
     else if(rowdif < 0) 
       rv = 1 ;
    else if(LocPos1.Z()>LocPos2.Z()) 
      rv = -1 ;
    else 
      rv = 1 ; 
     }

  else {
    if(PHOSMod1 < PHOSMod2 ) 
      rv = -1 ;
    else 
      rv = 1 ;
  }

  return rv ; 
}

//______________________________________________________________________________
void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
{
  //Execute action corresponding to one event
  //  This member function is called when a AliPHOSRecPoint is clicked with the locator
  //
  //  If Left button is clicked on AliPHOSRecPoint, the digits are switched on    
  //  and switched off when the mouse button is released.
  //

  //   static Int_t pxold, pyold;

   static TGraph *  DigitGraph = 0 ;

   if (!gPad->IsEditable()) return;

   TH2F * Histo = 0 ;
   TCanvas * HistoCanvas ; 
   
   switch (event) {
   
   case kButton1Down: {
     AliPHOSDigit * digit ;
     AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
     Int_t iDigit;
     Int_t relid[4] ;
     Float_t xi[fMulDigit] ;
     Float_t zi[fMulDigit] ;

     // create the histogram for the single cluster 
     // 1. gets histogram boundaries
     Float_t ximax = -999. ; 
     Float_t zimax = -999. ; 
     Float_t ximin = 999. ; 
     Float_t zimin = 999. ;
 
     for(iDigit=0; iDigit<fMulDigit; iDigit++) {
       digit = (AliPHOSDigit *) fDigitsList[iDigit];
       PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
       PHOSGeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
       if ( xi[iDigit] > ximax )
	 ximax = xi[iDigit] ; 
       if ( xi[iDigit] < ximin )
	 ximin = xi[iDigit] ; 
       if ( zi[iDigit] > zimax )
	 zimax = zi[iDigit] ; 
       if ( zi[iDigit] < zimin )
	 zimin = zi[iDigit] ;     
     }
     ximax += PHOSGeom->GetCrystalSize(0) / 2. ;
     zimax += PHOSGeom->GetCrystalSize(2) / 2. ;
     ximin -= PHOSGeom->GetCrystalSize(0) / 2. ;
     zimin -= PHOSGeom->GetCrystalSize(2) / 2. ;
     Int_t xdim = (int)( (ximax - ximin ) / PHOSGeom->GetCrystalSize(0) + 0.5  ) ; 
     Int_t zdim = (int)( (zimax - zimin ) / PHOSGeom->GetCrystalSize(2) + 0.5 ) ;
 
     // 2. gets the histogram title

     Text_t title[100] ; 
     sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
  
     if (!Histo) {
       delete Histo ; 
       Histo = 0 ; 
     }
     Histo = new TH2F("cluster3D", title,  xdim, ximin, ximax, zdim, zimin, zimax)  ;

     Float_t x, z ; 
     for(iDigit=0; iDigit<fMulDigit; iDigit++) {
       digit = (AliPHOSDigit *) fDigitsList[iDigit];
       PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
       PHOSGeom->RelPosInModule(relid, x, z);
       Histo->Fill(x, z, fEnergyList[iDigit] ) ;
     }

     if (!DigitGraph) {
       DigitGraph = new TGraph(fMulDigit,xi,zi);
       DigitGraph-> SetMarkerStyle(5) ; 
       DigitGraph-> SetMarkerSize(1.) ;
       DigitGraph-> SetMarkerColor(1) ;
       DigitGraph-> Paint("P") ;
     }

     Print() ;
     HistoCanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; 
     HistoCanvas->Draw() ; 
     Histo->Draw("lego1") ; 

     break;
   }

   case kButton1Up: 
     if (DigitGraph) {
       delete DigitGraph  ;
       DigitGraph = 0 ;
     }
     break;
  
   }
}

//____________________________________________________________________________
Float_t  AliPHOSEmcRecPoint::GetDispersion() 
{
  Float_t D    = 0 ;
  Float_t wtot = 0 ;

  TVector3 LocPos;
  GetLocalPosition(LocPos);
  Float_t x = LocPos.X() ;
  Float_t z = LocPos.Z() ;
  //  Int_t i = GetPHOSMod() ;

  AliPHOSDigit * digit ;
  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
  
  Int_t iDigit;
  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
    digit = (AliPHOSDigit *) fDigitsList[iDigit];
    Int_t relid[4] ;
    Float_t xi ;
    Float_t zi ;
    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
    PHOSGeom->RelPosInModule(relid, xi, zi);
    Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
    D += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; 
    wtot+=w ;
  }

  D /= wtot ;

  return TMath::Sqrt(D) ;
}

//____________________________________________________________________________
void  AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
{
  Float_t wtot = 0. ;
  Float_t x    = 0.;
  Float_t z    = 0.;
  Float_t Dxx  = 0.;
  Float_t Dzz  = 0.;
  Float_t Dxz  = 0.;

  AliPHOSDigit * digit ;
  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;
  Int_t iDigit;

  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
    digit = (AliPHOSDigit *) fDigitsList[iDigit];
    Int_t relid[4] ;
    Float_t xi ;
    Float_t zi ;
    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
    PHOSGeom->RelPosInModule(relid, xi, zi);
    Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
    Dxx  += w * xi * xi ;
    x    += w * xi ;
    Dzz  += w * zi * zi ;
    z    += w * zi ; 
    Dxz  += w * xi * zi ; 
    wtot += w ;
  }

  Dxx /= wtot ;
  x   /= wtot ;
  Dxx -= x * x ;
  Dzz /= wtot ;
  z   /= wtot ;
  Dzz -= z * z ;
  Dxz /= wtot ;
  Dxz -= x * z ;

  lambda[0] = TMath::Sqrt( 0.5 * (Dxx + Dzz) + TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
  lambda[1] = TMath::Sqrt( 0.5 * (Dxx + Dzz) - TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
}

//____________________________________________________________________________
Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
{
  Float_t menergy = 0. ;

  Int_t iDigit;

  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
 
    if(fEnergyList[iDigit] > menergy) 
      menergy = fEnergyList[iDigit] ;
  }
  return menergy ;
}

//____________________________________________________________________________
Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) 
{
  Int_t multipl   = 0 ;
  Int_t iDigit ;
  for(iDigit=0; iDigit<fMulDigit; iDigit++) {

    if(fEnergyList[iDigit] > H * fAmp) 
      multipl++ ;
  }
  return multipl ;
}

//____________________________________________________________________________
Int_t  AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t *  maxAt, Float_t * maxAtEnergy) 
{ 
  AliPHOSDigit * digit ;
  AliPHOSDigit * digitN ;
  

  Int_t iDigitN ;
  Int_t iDigit ;

  for(iDigit = 0; iDigit < fMulDigit; iDigit++){
    maxAt[iDigit] = fDigitsList[iDigit] ;
  }
  
  for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {   
    if(maxAt[iDigit] != -1) {
      digit = (AliPHOSDigit *) maxAt[iDigit] ;
         
      for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {	
	digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ; 
	
	if ( AreNeighbours(digit, digitN) ) {
	  if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {    
	    maxAt[iDigitN] = -1 ;
	    // but may be digit too is not local max ?
	    if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut) 
	      maxAt[iDigit] = -1 ;
	  }
	  else {
	    maxAt[iDigit] = -1 ;
	    // but may be digitN too is not local max ?
	    if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut) 
	      maxAt[iDigitN] = -1 ; 
	  } 
	} // if Areneighbours
      } // while digitN
    } // slot not empty
  } // while digit
  
  iDigitN = 0 ;
  for(iDigit = 0; iDigit < fMulDigit; iDigit++) { 
    if(maxAt[iDigit] != -1){
      maxAt[iDigitN] = maxAt[iDigit] ;
      maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
      iDigitN++ ; 
    }
  }
  return iDigitN ;
}

//____________________________________________________________________________
void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
{
  if( fLocPos.X() < 1000000.) { // already evaluated
   LPos = fLocPos ;
   return ;
  }

  Float_t wtot = 0. ;
 
  Int_t relid[4] ;

  Float_t x = 0. ;
  Float_t z = 0. ;
  
  AliPHOSDigit * digit ;

  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;

  Int_t iDigit;


  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
    digit = (AliPHOSDigit *) fDigitsList[iDigit];

    Float_t xi ;
    Float_t zi ;
    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
    PHOSGeom->RelPosInModule(relid, xi, zi);
    Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
    x    += xi * w ;
    z    += zi * w ;
    wtot += w ;

  }

  x /= wtot ;
  z /= wtot ;
  fLocPos.SetX(x)  ;
  fLocPos.SetY(0.) ;
  fLocPos.SetZ(z)  ;

  LPos = fLocPos ;
}

// //____________________________________________________________________________
// AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu) 
// {
//   int * DL = Clu.GetDigitsList() ; 
 
//  if(fDigitsList) 
//     delete fDigitsList  ;

//   AliPHOSDigit * digit ;
 
//   Int_t iDigit;

//   for(iDigit=0; iDigit<fMulDigit; iDigit++) {
//     digit = (AliPHOSDigit *) DL[iDigit];
//     AddDigit(*digit) ;
//   }

//   fAmp       = Clu.GetTotalEnergy() ;
//   fGeom      = Clu.GetGeom() ;
//   TVector3 LocPos;
//   Clu.GetLocalPosition(LocPos) ;
//   fLocPos    = LocPos;
//   fMulDigit       = Clu.GetMultiplicity() ;
//   fMaxDigit  = Clu.GetMaximumMultiplicity() ;
//   fPHOSMod   = Clu.GetPHOSMod() ;
//   fW0        = Clu.GetLogWeightCut() ; 
//   fDelta     = Clu.GetDelta() ;
//   fLocMaxCut = Clu.GetLocMaxCut() ;
  
//   delete DL ; 
 
//   return *this ;
// }

//____________________________________________________________________________
void AliPHOSEmcRecPoint::Print(Option_t * option) 
{
  cout << "AliPHOSEmcRecPoint: " << endl ;

  AliPHOSDigit * digit ; 
  Int_t iDigit;
  AliPHOSGeometry * PHOSGeom =  (AliPHOSGeometry *) fGeom ;

  Float_t xi ;
  Float_t zi ;
  Int_t relid[4] ; 
 
  for(iDigit=0; iDigit<fMulDigit; iDigit++) {
    digit = (AliPHOSDigit *) fDigitsList[iDigit];
    PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
    PHOSGeom->RelPosInModule(relid, xi, zi);
    cout << " Id = " << digit->GetId() ;  
    cout << "   module  = " << relid[0] ;  
    cout << "   x  = " << xi ;  
    cout << "   z  = " << zi ;  
    cout << "   Energy = " << fEnergyList[iDigit] << endl ;
  }
  cout << "       Multiplicity    = " << fMulDigit  << endl ;
  cout << "       Cluster Energy  = " << fAmp << endl ;
  
}

//______________________________________________________________________________
void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
{
   // Stream an object of class AliPHOSEmcRecPoint.

   if (R__b.IsReading()) {
      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
      AliPHOSRecPoint::Streamer(R__b);
      R__b >> fDelta;
      R__b >> fLocMaxCut;
      R__b.ReadArray(fEnergyList);
      R__b >> fW0;
   } else {
      R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
      AliPHOSRecPoint::Streamer(R__b);
      R__b << fDelta;
      R__b << fLocMaxCut;
      R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
      R__b << fW0;
   }
}
Commit	Line	Data
d15a28e7	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
	16	//_________________________________________________________________________
	17	// Rec Point in the PHOS EM calorimeter
	18	//*-- Author : Dmitri Peressounko RRC KI
	19	//////////////////////////////////////////////////////////////////////////////
	20
	21	// --- ROOT system ---
9f616d61	22	#include "TPad.h"
9f616d61	23	#include "TH2.h"
d15a28e7	24	#include "TMath.h"
9f616d61	25	#include "TCanvas.h"
d15a28e7	26
	27	// --- Standard library ---
	28
9f616d61	29	#include <iostream>
d15a28e7	30
	31	// --- AliRoot header files ---
	32
	33	#include "AliPHOSGeometry.h"
	34	#include "AliPHOSEmcRecPoint.h"
	35	#include "AliRun.h"
	36
	37	ClassImp(AliPHOSEmcRecPoint)
	38
	39	//____________________________________________________________________________
	40	AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut)
	41	: AliPHOSRecPoint()
	42	{
	43	// ctor
	44
	45	fMulDigit = 0 ;
	46	fAmp = 0. ;
	47	fEnergyList = new Float_t[fMaxDigit];
	48	AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
	49	fDelta = PHOSGeom->GetCrystalSize(0) ;
	50	fW0 = W0 ;
	51	fLocMaxCut = LocMaxCut ;
	52	fLocPos.SetX(1000000.) ; //Local position should be evaluated
	53	}
	54
9f616d61	55	//____________________________________________________________________________
	56	AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint()
	57	{
	58	// dtor
	59	}
d15a28e7	60
	61	//____________________________________________________________________________
	62	void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy)
	63	{
	64	// adds a digit to the digits list
	65	// and accumulates the total amplitude and the multiplicity
	66
	67	if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists
9f616d61	68	fMaxDigit*=2 ;
	69	int * tempo = new ( int[fMaxDigit] ) ;
	70	Float_t * tempoE = new ( Float_t[fMaxDigit] ) ;
	71
	72	Int_t index ;
d15a28e7	73	for ( index = 0 ; index < fMulDigit ; index++ ){
	74	tempo[index] = fDigitsList[index] ;
	75	tempoE[index] = fEnergyList[index] ;
	76	}
	77
9f616d61	78	delete [] fDigitsList ;
	79	fDigitsList = new ( int[fMaxDigit] ) ;
	80
	81	delete [] fEnergyList ;
	82	fEnergyList = new ( Float_t[fMaxDigit] ) ;
	83
	84	for ( index = 0 ; index < fMulDigit ; index++ ){
	85	fDigitsList[index] = tempo[index] ;
	86	fEnergyList[index] = tempoE[index] ;
	87	}
	88
	89	delete [] tempo ;
	90	delete [] tempoE ;
	91	} // if
d15a28e7	92
9f616d61	93	fDigitsList[fMulDigit] = (int) &digit ;
9f616d61	94	fEnergyList[fMulDigit++] = Energy ;
d15a28e7	95	fAmp += Energy ;
	96	}
	97
	98	//____________________________________________________________________________
	99	Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 )
	100	{
	101
	102	Bool_t aren = kFALSE ;
	103
	104	AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
	105	Int_t relid1[4] ;
	106	PHOSGeom->AbsToRelNumbering(digit1->GetId(), relid1) ;
	107
	108	Int_t relid2[4] ;
	109	PHOSGeom->AbsToRelNumbering(digit2->GetId(), relid2) ;
	110
	111	Int_t RowDiff = TMath::Abs( relid1[2] - relid2[2] ) ;
	112	Int_t ColDiff = TMath::Abs( relid1[3] - relid2[3] ) ;
	113
	114	if (( ColDiff<=1 ) && ( RowDiff <= 1 ) && (ColDiff+RowDiff > 0))
	115	aren = kTRUE ;
	116
	117	return aren ;
	118	}
	119
	120	//____________________________________________________________________________
	121	Int_t AliPHOSEmcRecPoint::Compare(TObject * obj)
	122	{
	123	Int_t rv ;
	124
	125	AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ;
	126
	127
	128	Int_t PHOSMod1 = this->GetPHOSMod() ;
	129	Int_t PHOSMod2 = clu->GetPHOSMod() ;
	130
	131	TVector3 LocPos1;
	132	this->GetLocalPosition(LocPos1) ;
	133	TVector3 LocPos2;
	134	clu->GetLocalPosition(LocPos2) ;
	135
	136	if(PHOSMod1 == PHOSMod2 ) {
	137	Int_t rowdif = (Int_t)TMath::Ceil(LocPos1.X()/fDelta)-(Int_t)TMath::Ceil(LocPos2.X()/fDelta) ;
	138	if (rowdif> 0)
	139	rv = -1 ;
	140	else if(rowdif < 0)
	141	rv = 1 ;
	142	else if(LocPos1.Z()>LocPos2.Z())
	143	rv = -1 ;
	144	else
	145	rv = 1 ;
	146	}
	147
	148	else {
	149	if(PHOSMod1 < PHOSMod2 )
	150	rv = -1 ;
	151	else
	152	rv = 1 ;
	153	}
	154
	155	return rv ;
	156	}
	157
9f616d61	158	//______________________________________________________________________________
	159	void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py)
	160	{
6ad0bfa0	161	//Execute action corresponding to one event
	162	// This member function is called when a AliPHOSRecPoint is clicked with the locator
	163	//
	164	// If Left button is clicked on AliPHOSRecPoint, the digits are switched on
	165	// and switched off when the mouse button is released.
	166	//
9f616d61	167
	168	// static Int_t pxold, pyold;
	169
	170	static TGraph * DigitGraph = 0 ;
	171
	172	if (!gPad->IsEditable()) return;
	173
	174	TH2F * Histo = 0 ;
	175	TCanvas * HistoCanvas ;
	176
	177	switch (event) {
	178
	179	case kButton1Down: {
	180	AliPHOSDigit * digit ;
	181	AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
	182	Int_t iDigit;
	183	Int_t relid[4] ;
	184	Float_t xi[fMulDigit] ;
	185	Float_t zi[fMulDigit] ;
	186
	187	// create the histogram for the single cluster
	188	// 1. gets histogram boundaries
	189	Float_t ximax = -999. ;
	190	Float_t zimax = -999. ;
	191	Float_t ximin = 999. ;
	192	Float_t zimin = 999. ;
	193
	194	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
	195	digit = (AliPHOSDigit *) fDigitsList[iDigit];
	196	PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
	197	PHOSGeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
	198	if ( xi[iDigit] > ximax )
	199	ximax = xi[iDigit] ;
	200	if ( xi[iDigit] < ximin )
	201	ximin = xi[iDigit] ;
	202	if ( zi[iDigit] > zimax )
	203	zimax = zi[iDigit] ;
	204	if ( zi[iDigit] < zimin )
	205	zimin = zi[iDigit] ;
	206	}
	207	ximax += PHOSGeom->GetCrystalSize(0) / 2. ;
	208	zimax += PHOSGeom->GetCrystalSize(2) / 2. ;
	209	ximin -= PHOSGeom->GetCrystalSize(0) / 2. ;
	210	zimin -= PHOSGeom->GetCrystalSize(2) / 2. ;
	211	Int_t xdim = (int)( (ximax - ximin ) / PHOSGeom->GetCrystalSize(0) + 0.5 ) ;
	212	Int_t zdim = (int)( (zimax - zimin ) / PHOSGeom->GetCrystalSize(2) + 0.5 ) ;
	213
	214	// 2. gets the histogram title
	215
	216	Text_t title[100] ;
	217	sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ;
	218
	219	if (!Histo) {
	220	delete Histo ;
	221	Histo = 0 ;
	222	}
	223	Histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ;
	224
	225	Float_t x, z ;
	226	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
	227	digit = (AliPHOSDigit *) fDigitsList[iDigit];
	228	PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
	229	PHOSGeom->RelPosInModule(relid, x, z);
	230	Histo->Fill(x, z, fEnergyList[iDigit] ) ;
231	}
232
233	if (!DigitGraph) {
234	DigitGraph = new TGraph(fMulDigit,xi,zi);
235	DigitGraph-> SetMarkerStyle(5) ;
236	DigitGraph-> SetMarkerSize(1.) ;
237	DigitGraph-> SetMarkerColor(1) ;
238	DigitGraph-> Paint("P") ;
239	}
240
241	Print() ;
242	HistoCanvas = new TCanvas("cluser", "a single cluster", 600, 500) ;
243	HistoCanvas->Draw() ;
244	Histo->Draw("lego1") ;
245
246	break;
247	}
248
249	case kButton1Up:
250	if (DigitGraph) {
251	delete DigitGraph ;
252	DigitGraph = 0 ;
253	}
254	break;
255
256	}
257	}
258
d15a28e7	259	//____________________________________________________________________________
	260	Float_t AliPHOSEmcRecPoint::GetDispersion()
	261	{
	262	Float_t D = 0 ;
	263	Float_t wtot = 0 ;
	264
	265	TVector3 LocPos;
	266	GetLocalPosition(LocPos);
	267	Float_t x = LocPos.X() ;
	268	Float_t z = LocPos.Z() ;
	269	// Int_t i = GetPHOSMod() ;
	270
	271	AliPHOSDigit * digit ;
	272	AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
	273
	274	Int_t iDigit;
	275	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
	276	digit = (AliPHOSDigit *) fDigitsList[iDigit];
	277	Int_t relid[4] ;
	278	Float_t xi ;
	279	Float_t zi ;
	280	PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
	281	PHOSGeom->RelPosInModule(relid, xi, zi);
	282	Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
	283	D += w((xi-x)(xi-x) + (zi-z)*(zi-z) ) ;
	284	wtot+=w ;
	285	}
	286
	287	D /= wtot ;
	288
	289	return TMath::Sqrt(D) ;
	290	}
	291
	292	//____________________________________________________________________________
	293	void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda)
	294	{
	295	Float_t wtot = 0. ;
	296	Float_t x = 0.;
	297	Float_t z = 0.;
	298	Float_t Dxx = 0.;
	299	Float_t Dzz = 0.;
	300	Float_t Dxz = 0.;
	301
	302	AliPHOSDigit * digit ;
	303	AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
	304	Int_t iDigit;
	305
	306	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
	307	digit = (AliPHOSDigit *) fDigitsList[iDigit];
	308	Int_t relid[4] ;
	309	Float_t xi ;
	310	Float_t zi ;
	311	PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
	312	PHOSGeom->RelPosInModule(relid, xi, zi);
	313	Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
	314	Dxx += w * xi * xi ;
	315	x += w * xi ;
	316	Dzz += w * zi * zi ;
	317	z += w * zi ;
	318	Dxz += w * xi * zi ;
	319	wtot += w ;
	320	}
	321
	322	Dxx /= wtot ;
323	x /= wtot ;
324	Dxx -= x * x ;
325	Dzz /= wtot ;
326	z /= wtot ;
327	Dzz -= z * z ;
328	Dxz /= wtot ;
329	Dxz -= x * z ;
330
331	lambda[0] = TMath::Sqrt( 0.5 * (Dxx + Dzz) + TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
332	lambda[1] = TMath::Sqrt( 0.5 * (Dxx + Dzz) - TMath::Sqrt( 0.25 * (Dxx - Dzz) * (Dxx - Dzz) + Dxz * Dxz ) ) ;
333	}
334
335	//____________________________________________________________________________
336	Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void)
337	{
338	Float_t menergy = 0. ;
339
340	Int_t iDigit;
341
342	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
343
344	if(fEnergyList[iDigit] > menergy)
345	menergy = fEnergyList[iDigit] ;
346	}
347	return menergy ;
348	}
349
350	//____________________________________________________________________________
351	Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H)
352	{
353	Int_t multipl = 0 ;
354	Int_t iDigit ;
355	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
356
357	if(fEnergyList[iDigit] > H * fAmp)
358	multipl++ ;
359	}
360	return multipl ;
361	}
362
363	//____________________________________________________________________________
6ad0bfa0	364	Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEnergy)
d15a28e7	365	{
	366	AliPHOSDigit * digit ;
	367	AliPHOSDigit * digitN ;
	368
	369
	370	Int_t iDigitN ;
	371	Int_t iDigit ;
	372
6ad0bfa0	373	for(iDigit = 0; iDigit < fMulDigit; iDigit++){
d15a28e7	374	maxAt[iDigit] = fDigitsList[iDigit] ;
	375	}
	376
6ad0bfa0	377	for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) {
9f616d61	378	if(maxAt[iDigit] != -1) {
9f616d61	379	digit = (AliPHOSDigit *) maxAt[iDigit] ;
d15a28e7	380
6ad0bfa0	381	for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) {
9f616d61	382	digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ;
d15a28e7	383
9f616d61	384	if ( AreNeighbours(digit, digitN) ) {
d15a28e7	385	if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) {
d15a28e7	386	maxAt[iDigitN] = -1 ;
6ad0bfa0	387	// but may be digit too is not local max ?
6ad0bfa0	388	if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut)
d15a28e7	389	maxAt[iDigit] = -1 ;
	390	}
	391	else {
	392	maxAt[iDigit] = -1 ;
6ad0bfa0	393	// but may be digitN too is not local max ?
6ad0bfa0	394	if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut)
d15a28e7	395	maxAt[iDigitN] = -1 ;
	396	}
	397	} // if Areneighbours
	398	} // while digitN
	399	} // slot not empty
	400	} // while digit
	401
	402	iDigitN = 0 ;
6ad0bfa0	403	for(iDigit = 0; iDigit < fMulDigit; iDigit++) {
d15a28e7	404	if(maxAt[iDigit] != -1){
d15a28e7	405	maxAt[iDigitN] = maxAt[iDigit] ;
9f616d61	406	maxAtEnergy[iDigitN] = fEnergyList[iDigit] ;
9f616d61	407	iDigitN++ ;
d15a28e7	408	}
	409	}
	410	return iDigitN ;
	411	}
	412
	413	//____________________________________________________________________________
	414	void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos)
	415	{
6ad0bfa0	416	if( fLocPos.X() < 1000000.) { // already evaluated
d15a28e7	417	LPos = fLocPos ;
	418	return ;
	419	}
	420
	421	Float_t wtot = 0. ;
	422
	423	Int_t relid[4] ;
	424
	425	Float_t x = 0. ;
	426	Float_t z = 0. ;
	427
	428	AliPHOSDigit * digit ;
	429
	430	AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
	431
	432	Int_t iDigit;
	433
	434
	435	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
	436	digit = (AliPHOSDigit *) fDigitsList[iDigit];
	437
	438	Float_t xi ;
	439	Float_t zi ;
	440	PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
	441	PHOSGeom->RelPosInModule(relid, xi, zi);
	442	Float_t w = TMath::Max( 0., fW0 + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ;
	443	x += xi * w ;
	444	z += zi * w ;
	445	wtot += w ;
	446
	447	}
	448
	449	x /= wtot ;
	450	z /= wtot ;
	451	fLocPos.SetX(x) ;
	452	fLocPos.SetY(0.) ;
	453	fLocPos.SetZ(z) ;
	454
	455	LPos = fLocPos ;
	456	}
	457
	458	// //____________________________________________________________________________
	459	// AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu)
	460	// {
	461	// int * DL = Clu.GetDigitsList() ;
	462
	463	// if(fDigitsList)
	464	// delete fDigitsList ;
	465
	466	// AliPHOSDigit * digit ;
	467
	468	// Int_t iDigit;
	469
	470	// for(iDigit=0; iDigit<fMulDigit; iDigit++) {
	471	// digit = (AliPHOSDigit *) DL[iDigit];
	472	// AddDigit(*digit) ;
	473	// }
	474
	475	// fAmp = Clu.GetTotalEnergy() ;
	476	// fGeom = Clu.GetGeom() ;
	477	// TVector3 LocPos;
	478	// Clu.GetLocalPosition(LocPos) ;
	479	// fLocPos = LocPos;
	480	// fMulDigit = Clu.GetMultiplicity() ;
481	// fMaxDigit = Clu.GetMaximumMultiplicity() ;
482	// fPHOSMod = Clu.GetPHOSMod() ;
483	// fW0 = Clu.GetLogWeightCut() ;
484	// fDelta = Clu.GetDelta() ;
485	// fLocMaxCut = Clu.GetLocMaxCut() ;
486
487	// delete DL ;
488
489	// return *this ;
490	// }
491
492	//____________________________________________________________________________
493	void AliPHOSEmcRecPoint::Print(Option_t * option)
494	{
495	cout << "AliPHOSEmcRecPoint: " << endl ;
496
497	AliPHOSDigit * digit ;
498	Int_t iDigit;
9f616d61	499	AliPHOSGeometry * PHOSGeom = (AliPHOSGeometry *) fGeom ;
d15a28e7	500
9f616d61	501	Float_t xi ;
	502	Float_t zi ;
	503	Int_t relid[4] ;
	504
d15a28e7	505	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
d15a28e7	506	digit = (AliPHOSDigit *) fDigitsList[iDigit];
9f616d61	507	PHOSGeom->AbsToRelNumbering(digit->GetId(), relid) ;
	508	PHOSGeom->RelPosInModule(relid, xi, zi);
	509	cout << " Id = " << digit->GetId() ;
	510	cout << " module = " << relid[0] ;
	511	cout << " x = " << xi ;
	512	cout << " z = " << zi ;
	513	cout << " Energy = " << fEnergyList[iDigit] << endl ;
d15a28e7	514	}
	515	cout << " Multiplicity = " << fMulDigit << endl ;
	516	cout << " Cluster Energy = " << fAmp << endl ;
	517
	518	}
	519
	520	//______________________________________________________________________________
	521	void AliPHOSEmcRecPoint::Streamer(TBuffer &R__b)
	522	{
	523	// Stream an object of class AliPHOSEmcRecPoint.
	524
	525	if (R__b.IsReading()) {
	526	Version_t R__v = R__b.ReadVersion(); if (R__v) { }
	527	AliPHOSRecPoint::Streamer(R__b);
	528	R__b >> fDelta;
	529	R__b >> fLocMaxCut;
	530	R__b.ReadArray(fEnergyList);
	531	R__b >> fW0;
	532	} else {
	533	R__b.WriteVersion(AliPHOSEmcRecPoint::IsA());
	534	AliPHOSRecPoint::Streamer(R__b);
	535	R__b << fDelta;
	536	R__b << fLocMaxCut;
	537	R__b.WriteArray(fEnergyList, GetMaximumDigitMultiplicity() );
	538	R__b << fW0;
	539	}
	540	}
	541
	542