[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];
92862013	48	AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
92862013	49	fDelta = phosgeom->GetCrystalSize(0) ;
d15a28e7	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
92862013	104	AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7	105	Int_t relid1[4] ;
92862013	106	phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ;
d15a28e7	107
d15a28e7	108	Int_t relid2[4] ;
92862013	109	phosgeom->AbsToRelNumbering(digit2->GetId(), relid2) ;
d15a28e7	110
92862013	111	Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
92862013	112	Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
d15a28e7	113
92862013	114	if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0))
d15a28e7	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
92862013	128	Int_t phosmod1 = this->GetPHOSMod() ;
92862013	129	Int_t phosmod2 = clu->GetPHOSMod() ;
d15a28e7	130
92862013	131	TVector3 locpos1;
	132	this->GetLocalPosition(locpos1) ;
	133	TVector3 locpos2;
	134	clu->GetLocalPosition(locpos2) ;
d15a28e7	135
92862013	136	if(phosmod1 == phosmod2 ) {
92862013	137	Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/fDelta)-(Int_t)TMath::Ceil(locpos2.X()/fDelta) ;
d15a28e7	138	if (rowdif> 0)
	139	rv = -1 ;
	140	else if(rowdif < 0)
	141	rv = 1 ;
92862013	142	else if(locpos1.Z()>locpos2.Z())
d15a28e7	143	rv = -1 ;
	144	else
	145	rv = 1 ;
	146	}
	147
	148	else {
92862013	149	if(phosmod1 < phosmod2 )
d15a28e7	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
92862013	170	static TGraph * digitgraph = 0 ;
9f616d61	171
	172	if (!gPad->IsEditable()) return;
	173
92862013	174	TH2F * histo = 0 ;
92862013	175	TCanvas * histocanvas ;
9f616d61	176
	177	switch (event) {
	178
	179	case kButton1Down: {
	180	AliPHOSDigit * digit ;
92862013	181	AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
9f616d61	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];
92862013	196	phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
92862013	197	phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]);
9f616d61	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	}
92862013	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 ) ;
9f616d61	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
92862013	219	if (!histo) {
	220	delete histo ;
	221	histo = 0 ;
9f616d61	222	}
92862013	223	histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ;
9f616d61	224
	225	Float_t x, z ;
	226	for(iDigit=0; iDigit<fMulDigit; iDigit++) {
	227	digit = (AliPHOSDigit *) fDigitsList[iDigit];
92862013	228	phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
	229	phosgeom->RelPosInModule(relid, x, z);
	230	histo->Fill(x, z, fEnergyList[iDigit] ) ;
9f616d61	231	}
9f616d61	232
92862013	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") ;
9f616d61	239	}
	240
	241	Print() ;
92862013	242	histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ;
	243	histocanvas->Draw() ;
	244	histo->Draw("lego1") ;
9f616d61	245
	246	break;
	247	}
	248
	249	case kButton1Up:
92862013	250	if (digitgraph) {
	251	delete digitgraph ;
	252	digitgraph = 0 ;
9f616d61	253	}
	254	break;
	255
	256	}
	257	}
	258
d15a28e7	259	//____________________________________________________________________________
	260	Float_t AliPHOSEmcRecPoint::GetDispersion()
	261	{
92862013	262	Float_t d = 0 ;
d15a28e7	263	Float_t wtot = 0 ;
d15a28e7	264
92862013	265	TVector3 locpos;
	266	GetLocalPosition(locpos);
	267	Float_t x = locpos.X() ;
	268	Float_t z = locpos.Z() ;
d15a28e7	269	// Int_t i = GetPHOSMod() ;
	270
	271	AliPHOSDigit * digit ;
92862013	272	AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7	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 ;
92862013	280	phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
92862013	281	phosgeom->RelPosInModule(relid, xi, zi);
d15a28e7	282	Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
92862013	283	d += w((xi-x)(xi-x) + (zi-z)*(zi-z) ) ;
d15a28e7	284	wtot+=w ;
	285	}
	286
92862013	287	d /= wtot ;
d15a28e7	288
92862013	289	return TMath::Sqrt(d) ;
d15a28e7	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.;
92862013	298	Float_t dxx = 0.;
	299	Float_t dzz = 0.;
	300	Float_t dxz = 0.;
d15a28e7	301
d15a28e7	302	AliPHOSDigit * digit ;
92862013	303	AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7	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 ;
92862013	311	phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
92862013	312	phosgeom->RelPosInModule(relid, xi, zi);
d15a28e7	313	Float_t w = TMath::Max(0.,fW0+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ;
92862013	314	dxx += w * xi * xi ;
d15a28e7	315	x += w * xi ;
92862013	316	dzz += w * zi * zi ;
d15a28e7	317	z += w * zi ;
92862013	318	dxz += w * xi * zi ;
d15a28e7	319	wtot += w ;
	320	}
	321
92862013	322	dxx /= wtot ;
d15a28e7	323	x /= wtot ;
92862013	324	dxx -= x * x ;
92862013	325	dzz /= wtot ;
d15a28e7	326	z /= wtot ;
92862013	327	dzz -= z * z ;
	328	dxz /= wtot ;
	329	dxz -= x * z ;
d15a28e7	330
92862013	331	lambda[0] = TMath::Sqrt( 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
92862013	332	lambda[1] = TMath::Sqrt( 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ) ;
d15a28e7	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
92862013	430	AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ;
d15a28e7	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 ;
92862013	440	phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
92862013	441	phosgeom->RelPosInModule(relid, xi, zi);
d15a28e7	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	// {
92862013	461	// int * dl = Clu.GetDigitsList() ;
d15a28e7	462
	463	// if(fDigitsList)
	464	// delete fDigitsList ;
	465
	466	// AliPHOSDigit * digit ;
	467
	468	// Int_t iDigit;
	469
	470	// for(iDigit=0; iDigit<fMulDigit; iDigit++) {
92862013	471	// digit = (AliPHOSDigit *) dl[iDigit];
d15a28e7	472	// AddDigit(*digit) ;
	473	// }
	474
	475	// fAmp = Clu.GetTotalEnergy() ;
	476	// fGeom = Clu.GetGeom() ;
92862013	477	// TVector3 locpos;
	478	// Clu.GetLocalPosition(locpos) ;
	479	// fLocPos = locpos;
d15a28e7	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
92862013	487	// delete dl ;
d15a28e7	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;
92862013	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];
92862013	507	phosgeom->AbsToRelNumbering(digit->GetId(), relid) ;
92862013	508	phosgeom->RelPosInModule(relid, xi, zi);
9f616d61	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