X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSEmcRecPoint.cxx;h=468f2a4968ad28b0bac1db16b46a6a1ba6084a31;hb=9bd51fae56c47e270d0d24cc4fb5281f94c86312;hp=df8047056ea31f32670eb6eab24a250172bcd42f;hpb=9286201303faf22df30f835f6aa39b6517513db8;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSEmcRecPoint.cxx b/PHOS/AliPHOSEmcRecPoint.cxx index df8047056ea..468f2a4968a 100644 --- a/PHOS/AliPHOSEmcRecPoint.cxx +++ b/PHOS/AliPHOSEmcRecPoint.cxx @@ -13,10 +13,14 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ +/* $Id$ */ + //_________________________________________________________________________ -// Rec Point in the PHOS EM calorimeter -//*-- Author : Dmitri Peressounko RRC KI -////////////////////////////////////////////////////////////////////////////// +// RecPoint implementation for PHOS-EMC +// An EmcRecPoint is a cluster of digits +//*-- +//*-- Author: Dmitri Peressounko (RRC KI & SUBATECH) + // --- ROOT system --- #include "TPad.h" @@ -26,57 +30,75 @@ // --- Standard library --- -#include - // --- AliRoot header files --- + #include "AliGenerator.h" #include "AliPHOSGeometry.h" #include "AliPHOSEmcRecPoint.h" #include "AliRun.h" +#include "AliPHOSGetter.h" ClassImp(AliPHOSEmcRecPoint) //____________________________________________________________________________ -AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut) - : AliPHOSRecPoint() +AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : AliPHOSRecPoint() { // ctor fMulDigit = 0 ; fAmp = 0. ; - fEnergyList = new Float_t[fMaxDigit]; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; - fDelta = phosgeom->GetCrystalSize(0) ; - fW0 = W0 ; - fLocMaxCut = LocMaxCut ; + fCoreEnergy = 0 ; + fEnergyList = 0 ; + fTime = -1. ; fLocPos.SetX(1000000.) ; //Local position should be evaluated + } //____________________________________________________________________________ -AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() +AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const char * opt) : AliPHOSRecPoint(opt) { - // dtor + // ctor + + fMulDigit = 0 ; + fAmp = 0. ; + fCoreEnergy = 0 ; + fEnergyList = 0 ; + fTime = -1. ; + fLocPos.SetX(1000000.) ; //Local position should be evaluated + +} + +//____________________________________________________________________________ +AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() +{ + // dtor + + if ( fEnergyList ) + delete[] fEnergyList ; } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy) +void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy) { - // adds a digit to the digits list + // Adds a digit to the RecPoint // and accumulates the total amplitude and the multiplicity + if(fEnergyList == 0) + fEnergyList = new Float_t[fMaxDigit]; + if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists fMaxDigit*=2 ; - int * tempo = new ( int[fMaxDigit] ) ; + Int_t * tempo = new ( Int_t[fMaxDigit] ) ; Float_t * tempoE = new ( Float_t[fMaxDigit] ) ; Int_t index ; for ( index = 0 ; index < fMulDigit ; index++ ){ - tempo[index] = fDigitsList[index] ; + tempo[index] = fDigitsList[index] ; tempoE[index] = fEnergyList[index] ; } delete [] fDigitsList ; - fDigitsList = new ( int[fMaxDigit] ) ; + fDigitsList = new ( Int_t[fMaxDigit] ) ; delete [] fEnergyList ; fEnergyList = new ( Float_t[fMaxDigit] ) ; @@ -90,18 +112,24 @@ void AliPHOSEmcRecPoint::AddDigit(AliDigitNew & digit, Float_t Energy) delete [] tempoE ; } // if - fDigitsList[fMulDigit] = (int) &digit ; - fEnergyList[fMulDigit++] = Energy ; + fDigitsList[fMulDigit] = digit.GetIndexInList() ; + fEnergyList[fMulDigit] = Energy ; + fMulDigit++ ; fAmp += Energy ; + + EvalPHOSMod(&digit) ; } //____________________________________________________________________________ -Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) +Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) const { + // Tells if (true) or not (false) two digits are neighbors Bool_t aren = kFALSE ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; + AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; + AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + Int_t relid1[4] ; phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; @@ -118,27 +146,32 @@ Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * d } //____________________________________________________________________________ -Int_t AliPHOSEmcRecPoint::Compare(TObject * obj) +Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const { + // Compares two RecPoints according to their position in the PHOS modules + + Float_t delta = 1 ; //Width of "Sorting row". If you changibg this + //value (what is senseless) change as vell delta in + //AliPHOSTrackSegmentMakerv* and other RecPoints... Int_t rv ; AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; - Int_t phosmod1 = this->GetPHOSMod() ; + Int_t phosmod1 = GetPHOSMod() ; Int_t phosmod2 = clu->GetPHOSMod() ; TVector3 locpos1; - this->GetLocalPosition(locpos1) ; + 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) ; + Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ; if (rowdif> 0) - rv = -1 ; + rv = 1 ; else if(rowdif < 0) - rv = 1 ; + rv = -1 ; else if(locpos1.Z()>locpos2.Z()) rv = -1 ; else @@ -154,163 +187,255 @@ Int_t AliPHOSEmcRecPoint::Compare(TObject * obj) return rv ; } - //______________________________________________________________________________ -void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) +void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) const { - //Execute action corresponding to one event + + // 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. - // + + + AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; + if(!gime) return ; + AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + + static TGraph * digitgraph = 0 ; + + if (!gPad->IsEditable()) return; + + TH2F * histo = 0 ; + TCanvas * histocanvas ; - // static Int_t pxold, pyold; + const TClonesArray * digits = gime->Digits() ; + + switch (event) { + + case kButton1Down: { + AliPHOSDigit * digit ; + Int_t iDigit; + Int_t relid[4] ; + + const Int_t kMulDigit = AliPHOSEmcRecPoint::GetDigitsMultiplicity() ; + Float_t * xi = new Float_t[kMulDigit] ; + Float_t * zi = new Float_t[kMulDigit] ; + + // 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; iDigitAt(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; iDigitAt(fDigitsList[iDigit]) ; + phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; + phosgeom->RelPosInModule(relid, x, z); + histo->Fill(x, z, fEnergyList[iDigit] ) ; + } + + if (!digitgraph) { + digitgraph = new TGraph(kMulDigit,xi,zi); + digitgraph-> SetMarkerStyle(5) ; + digitgraph-> SetMarkerSize(1.) ; + digitgraph-> SetMarkerColor(1) ; + digitgraph-> Paint("P") ; + } + + // Print() ; + histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ; + histocanvas->Draw() ; + histo->Draw("lego1") ; + + delete[] xi ; + delete[] zi ; + + break; + } + + case kButton1Up: + if (digitgraph) { + delete digitgraph ; + digitgraph = 0 ; + } + break; + + } +} - static TGraph * digitgraph = 0 ; +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits) +{ + // Calculates the dispersion of the shower at the origine of the RecPoint - if (!gPad->IsEditable()) return; + Float_t d = 0. ; + Float_t wtot = 0. ; - 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; iDigitAbsToRelNumbering(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 ) ; + Float_t x = 0.; + Float_t z = 0.; + + AliPHOSDigit * digit ; - // 2. gets the histogram title + AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; + AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + - Text_t title[100] ; - sprintf(title,"Energy=%1.2f GeV ; Digits ; %d ", GetEnergy(), GetDigitsMultiplicity()) ; + // Calculates the center of gravity in the local PHOS-module coordinates - if (!histo) { - delete histo ; - histo = 0 ; - } - histo = new TH2F("cluster3D", title, xdim, ximin, ximax, zdim, zimin, zimax) ; - - Float_t x, z ; - for(iDigit=0; iDigitAbsToRelNumbering(digit->GetId(), relid) ; - phosgeom->RelPosInModule(relid, x, z); - histo->Fill(x, z, fEnergyList[iDigit] ) ; - } + Int_t iDigit; - if (!digitgraph) { - digitgraph = new TGraph(fMulDigit,xi,zi); - digitgraph-> SetMarkerStyle(5) ; - digitgraph-> SetMarkerSize(1.) ; - digitgraph-> SetMarkerColor(1) ; - digitgraph-> Paint("P") ; - } + for(iDigit=0; iDigitAt(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., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; + x += xi * w ; + z += zi * w ; + wtot += w ; + } + x /= wtot ; + z /= wtot ; - Print() ; - histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; - histocanvas->Draw() ; - histo->Draw("lego1") ; - break; - } +// Calculates the dispersion in coordinates + wtot = 0.; + for(iDigit=0; iDigit < fMulDigit; iDigit++) { + digit = (AliPHOSDigit *) digits->At(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.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; + wtot+=w ; - case kButton1Up: - if (digitgraph) { - delete digitgraph ; - digitgraph = 0 ; - } - break; - } -} + } + -//____________________________________________________________________________ -Float_t AliPHOSEmcRecPoint::GetDispersion() + d /= wtot ; + + fDispersion = TMath::Sqrt(d) ; + +} +//______________________________________________________________________________ +void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) { - Float_t d = 0 ; - Float_t wtot = 0 ; + // This function calculates energy in the core, + // i.e. within a radius rad = 3cm around the center. Beyond this radius + // in accordance with shower profile the energy deposition + // should be less than 2% - TVector3 locpos; - GetLocalPosition(locpos); - Float_t x = locpos.X() ; - Float_t z = locpos.Z() ; - // Int_t i = GetPHOSMod() ; + Float_t coreRadius = 3 ; + + Float_t x = 0 ; + Float_t z = 0 ; AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; - + + AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; + AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + Int_t iDigit; + +// Calculates the center of gravity in the local PHOS-module coordinates + Float_t wtot = 0; for(iDigit=0; iDigitAt(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 ; + Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; + x += xi * w ; + z += zi * w ; + wtot += w ; } + x /= wtot ; + z /= wtot ; - d /= wtot ; - return TMath::Sqrt(d) ; + for(iDigit=0; iDigit < fMulDigit; iDigit++) { + digit = (AliPHOSDigit *) ( digits->At(fDigitsList[iDigit]) ) ; + Int_t relid[4] ; + Float_t xi ; + Float_t zi ; + phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; + phosgeom->RelPosInModule(relid, xi, zi); + Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(zi-z)*(zi-z)) ; + if(distance < coreRadius) + fCoreEnergy += fEnergyList[iDigit] ; + } + } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda) +void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) { - Float_t wtot = 0. ; - Float_t x = 0.; - Float_t z = 0.; - Float_t dxx = 0.; - Float_t dzz = 0.; - Float_t dxz = 0.; + // Calculates the axis of the shower ellipsoid + + Double_t wtot = 0. ; + Double_t x = 0.; + Double_t z = 0.; + Double_t dxx = 0.; + Double_t dzz = 0.; + Double_t dxz = 0.; AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; + + AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; + AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + Int_t iDigit; + for(iDigit=0; iDigitAt(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 ) ) ; + Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; dxx += w * xi * xi ; x += w * xi ; dzz += w * zi * zi ; @@ -318,7 +443,6 @@ void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda) dxz += w * xi * zi ; wtot += w ; } - dxx /= wtot ; x /= wtot ; dxx -= x * x ; @@ -328,13 +452,113 @@ void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda) 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 ) ) ; +// //Apply correction due to non-perpendicular incidence +// Double_t CosX ; +// Double_t CosZ ; +// AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; +// AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + // Double_t DistanceToIP= (Double_t ) phosgeom->GetIPtoCrystalSurface() ; + +// CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ; +// CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ; + +// dxx = dxx/(CosX*CosX) ; +// dzz = dzz/(CosZ*CosZ) ; +// dxz = dxz/(CosX*CosZ) ; + + + fLambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; + if(fLambda[0] > 0) + fLambda[0] = TMath::Sqrt(fLambda[0]) ; + + fLambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; + if(fLambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda. + fLambda[1] = TMath::Sqrt(fLambda[1]) ; + else + fLambda[1]= 0. ; } //____________________________________________________________________________ -Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) +void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits ) +{ + // Evaluates all shower parameters + + AliPHOSRecPoint::EvalAll(logWeight,digits) ; + EvalLocalPosition(logWeight, digits) ; + EvalElipsAxis(logWeight, digits) ; + EvalDispersion(logWeight, digits) ; + EvalCoreEnergy(logWeight, digits); + EvalTime(digits) ; +} +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits) { + // Calculates the center of gravity in the local PHOS-module coordinates + Float_t wtot = 0. ; + + Int_t relid[4] ; + + Float_t x = 0. ; + Float_t z = 0. ; + + AliPHOSDigit * digit ; + + AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; + AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + + Int_t iDigit; + + for(iDigit=0; iDigitAt(fDigitsList[iDigit]) ; + + Float_t xi ; + Float_t zi ; + phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; + phosgeom->RelPosInModule(relid, xi, zi); + Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; + x += xi * w ; + z += zi * w ; + wtot += w ; + + } + + x /= wtot ; + z /= wtot ; + + // Correction for the depth of the shower starting point (TDR p 127) + Float_t para = 0.925 ; + Float_t parb = 6.52 ; + + Float_t xo,yo,zo ; //Coordinates of the origin + gAlice->Generator()->GetOrigin(xo,yo,zo) ; + + Float_t phi = phosgeom->GetPHOSAngle(relid[0]) ; + + //Transform to the local ref.frame + Float_t xoL,yoL ; + xoL = xo*TMath::Cos(phi)-yo*TMath::Sin(phi) ; + yoL = xo*TMath::Sin(phi)+yo*TMath::Cos(phi) ; + + Float_t radius = phosgeom->GetIPtoCrystalSurface()-yoL; + + Float_t incidencephi = TMath::ATan((x-xoL ) / radius) ; + Float_t incidencetheta = TMath::ATan((z-zo) / radius) ; + + Float_t depthx = ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencephi) ; + Float_t depthz = ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencetheta) ; + + fLocPos.SetX(x - depthx) ; + fLocPos.SetY(0.) ; + fLocPos.SetZ(z - depthz) ; + + fLocPosM = 0 ; +} + +//____________________________________________________________________________ +Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) const +{ + // Finds the maximum energy in the cluster + Float_t menergy = 0. ; Int_t iDigit; @@ -348,8 +572,10 @@ Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) } //____________________________________________________________________________ -Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) +Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) const { + // Calculates the multiplicity of digits with energy larger than H*energy + Int_t multipl = 0 ; Int_t iDigit ; for(iDigit=0; iDigitAt(fDigitsList[iDigit]) ; + for(iDigit = 0 ; iDigit < fMulDigit; iDigit++) { - if(maxAt[iDigit] != -1) { - digit = (AliPHOSDigit *) maxAt[iDigit] ; - + if(maxAt[iDigit]) { + digit = maxAt[iDigit] ; + for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) { - digitN = (AliPHOSDigit *) fDigitsList[iDigitN] ; + digitN = (AliPHOSDigit *) digits->At(fDigitsList[iDigitN]) ; if ( AreNeighbours(digit, digitN) ) { if (fEnergyList[iDigit] > fEnergyList[iDigitN] ) { - maxAt[iDigitN] = -1 ; + maxAt[iDigitN] = 0 ; // but may be digit too is not local max ? - if(fEnergyList[iDigit] < fEnergyList[iDigitN] + fLocMaxCut) - maxAt[iDigit] = -1 ; + if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut) + maxAt[iDigit] = 0 ; } else { - maxAt[iDigit] = -1 ; + maxAt[iDigit] = 0 ; // but may be digitN too is not local max ? - if(fEnergyList[iDigit] > fEnergyList[iDigitN] - fLocMaxCut) - maxAt[iDigitN] = -1 ; + if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut) + maxAt[iDigitN] = 0 ; } } // if Areneighbours } // while digitN @@ -401,7 +631,7 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEn iDigitN = 0 ; for(iDigit = 0; iDigit < fMulDigit; iDigit++) { - if(maxAt[iDigit] != -1){ + if(maxAt[iDigit]){ maxAt[iDigitN] = maxAt[iDigit] ; maxAtEnergy[iDigitN] = fEnergyList[iDigit] ; iDigitN++ ; @@ -409,134 +639,79 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEn } return iDigitN ; } - //____________________________________________________________________________ -void AliPHOSEmcRecPoint::GetLocalPosition(TVector3 &LPos) -{ - if( fLocPos.X() < 1000000.) { // already evaluated - LPos = fLocPos ; - return ; +void AliPHOSEmcRecPoint::EvalTime(TClonesArray * digits){ + + Float_t maxE = 0; + Int_t maxAt = 0; + for(Int_t idig=0; idig < fMulDigit; idig++){ + if(fEnergyList[idig] > maxE){ + maxE = fEnergyList[idig] ; + maxAt = idig; + } } - - Float_t wtot = 0. ; - - Int_t relid[4] ; - - Float_t x = 0. ; - Float_t z = 0. ; + fTime = ((AliPHOSDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ; - AliPHOSDigit * digit ; - - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; - - Int_t iDigit; - - - for(iDigit=0; iDigitAbsToRelNumbering(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 ; - +} +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::Purify(Float_t threshold){ + //Removes digits below threshold + + Int_t * tempo = new ( Int_t[fMaxDigit] ) ; + Float_t * tempoE = new ( Float_t[fMaxDigit] ) ; + + Int_t mult = 0 ; + for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){ + if(fEnergyList[iDigit] > threshold){ + tempo[mult] = fDigitsList[iDigit] ; + tempoE[mult] = fEnergyList[iDigit] ; + mult++ ; + } } + + fMulDigit = mult ; + delete [] fDigitsList ; + delete [] fEnergyList ; + fDigitsList = new (Int_t[fMulDigit]) ; + fEnergyList = new ( Float_t[fMulDigit]) ; + + for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){ + fDigitsList[iDigit] = tempo[iDigit]; + fEnergyList[iDigit] = tempoE[iDigit] ; + } + + delete [] tempo ; + delete [] tempoE ; - 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; iDigitAbsToRelNumbering(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 ; + for(iDigit=0; iDigit> 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; - } -} - - + +