X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=PHOS%2FAliPHOSEmcRecPoint.cxx;h=2c9625246a11da7a1ccf806bcc7f7c999aeb9133;hp=8d7d035613ce50affbdc610755020e620aacf564;hb=0488b68a2579fab9091b4daba6be38aef4002543;hpb=fc7e2f436398e9b9b8875125b25f6cb4f7e97591 diff --git a/PHOS/AliPHOSEmcRecPoint.cxx b/PHOS/AliPHOSEmcRecPoint.cxx index 8d7d035613c..2c9625246a1 100644 --- a/PHOS/AliPHOSEmcRecPoint.cxx +++ b/PHOS/AliPHOSEmcRecPoint.cxx @@ -15,6 +15,29 @@ /* $Id$ */ +/* History of cvs commits: + * + * $Log$ + * Revision 1.57 2007/04/05 10:18:58 policheh + * Introduced distance to nearest bad crystal. + * + * Revision 1.56 2007/03/06 06:47:28 kharlov + * DP:Possibility to use actual vertex position added + * + * Revision 1.55 2007/01/19 20:31:19 kharlov + * Improved formatting for Print() + * + * Revision 1.54 2006/08/28 10:01:56 kharlov + * Effective C++ warnings fixed (Timur Pocheptsov) + * + * Revision 1.53 2005/12/20 14:28:47 hristov + * Additional protection + * + * Revision 1.52 2005/05/28 14:19:04 schutz + * Compilation warnings fixed by T.P. + * + */ + //_________________________________________________________________________ // RecPoint implementation for PHOS-EMC // An EmcRecPoint is a cluster of digits @@ -31,6 +54,7 @@ // --- Standard library --- // --- AliRoot header files --- +#include "AliLog.h" #include "AliPHOSLoader.h" #include "AliGenerator.h" #include "AliPHOSGeometry.h" @@ -40,58 +64,54 @@ ClassImp(AliPHOSEmcRecPoint) //____________________________________________________________________________ -AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : AliPHOSRecPoint() +AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : + AliPHOSRecPoint(), + fCoreEnergy(0.), fDispersion(0.), + fEnergyList(0), fTime(-1.), fNExMax(0), + fM2x(0.), fM2z(0.), fM3x(0.), fM4z(0.), + fPhixe(0.), fDistToBadCrystal(-1),fDebug(0) { // ctor - fMulDigit = 0 ; fAmp = 0. ; - fCoreEnergy = 0 ; - fEnergyList = 0 ; - fNExMax = 0 ; //Not unfolded yet - fTime = -1. ; fLocPos.SetX(1000000.) ; //Local position should be evaluated - fDebug=0; - } //____________________________________________________________________________ -AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const char * opt) : AliPHOSRecPoint(opt) +AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const char * opt) : + AliPHOSRecPoint(opt), + fCoreEnergy(0.), fDispersion(0.), + fEnergyList(0), fTime(-1.), fNExMax(0), + fM2x(0.), fM2z(0.), fM3x(0.), fM4z(0.), + fPhixe(0.), fDistToBadCrystal(-1), fDebug(0) { // ctor - fMulDigit = 0 ; fAmp = 0. ; - fNExMax = 0 ; //Not unfolded yet - fCoreEnergy = 0 ; - fEnergyList = 0 ; - fTime = -1. ; fLocPos.SetX(1000000.) ; //Local position should be evaluated - fDebug=0; - } //____________________________________________________________________________ -AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const AliPHOSEmcRecPoint & rp) : AliPHOSRecPoint(rp) +AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const AliPHOSEmcRecPoint & rp) : + AliPHOSRecPoint(rp), + fCoreEnergy(rp.fCoreEnergy), fDispersion(rp.fDispersion), + fEnergyList(0), fTime(rp.fTime), fNExMax(rp.fNExMax), + fM2x(rp.fM2x), fM2z(rp.fM2z), fM3x(rp.fM3x), fM4z(rp.fM4z), + fPhixe(rp.fPhixe), fDistToBadCrystal(rp.fDistToBadCrystal), fDebug(rp.fDebug) { // cpy ctor - fMulDigit = rp.fMulDigit ; fAmp = rp.fAmp ; - fCoreEnergy = rp.fCoreEnergy ; fEnergyList = new Float_t[rp.fMulDigit] ; Int_t index ; for(index = 0 ; index < fMulDigit ; index++) fEnergyList[index] = rp.fEnergyList[index] ; - fNExMax = rp.fNExMax ; - fTime = rp.fTime ; } //____________________________________________________________________________ AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() { // dtor - if ( fEnergyList ) delete[] fEnergyList ; } @@ -107,8 +127,8 @@ void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy) if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists fMaxDigit*=2 ; - Int_t * tempo = new ( Int_t[fMaxDigit] ) ; - Float_t * tempoE = new ( Float_t[fMaxDigit] ) ; + Int_t * tempo = new Int_t[fMaxDigit]; + Float_t * tempoE = new Float_t[fMaxDigit]; Int_t index ; for ( index = 0 ; index < fMulDigit ; index++ ){ @@ -117,10 +137,10 @@ void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy) } delete [] fDigitsList ; - fDigitsList = new ( Int_t[fMaxDigit] ) ; + fDigitsList = new Int_t[fMaxDigit]; delete [] fEnergyList ; - fEnergyList = new ( Float_t[fMaxDigit] ) ; + fEnergyList = new Float_t[fMaxDigit]; for ( index = 0 ; index < fMulDigit ; index++ ){ fDigitsList[index] = tempo[index] ; @@ -146,7 +166,7 @@ Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * d Bool_t aren = kFALSE ; - AliPHOSGeometry * phosgeom = AliPHOSLoader::GetPHOSGeometry(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; Int_t relid1[4] ; phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; @@ -167,8 +187,8 @@ Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * d 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 + + const 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 ; @@ -206,7 +226,7 @@ Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const return rv ; } //______________________________________________________________________________ -void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) const +void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) /*const*/ { // Execute action corresponding to one event @@ -216,7 +236,7 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) const // and switched off when the mouse button is released. - AliPHOSGeometry * phosgeom = AliPHOSLoader::GetPHOSGeometry(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); static TGraph * digitgraph = 0 ; @@ -227,16 +247,16 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) const //try to get run loader from default event folder - AliRunLoader* rn = AliRunLoader::GetRunLoader(AliConfig::fgkDefaultEventFolderName); + AliRunLoader* rn = AliRunLoader::GetRunLoader(AliConfig::GetDefaultEventFolderName()); if (rn == 0x0) { - Error("ExecuteEvent","Can not find Run Loader in Default Event Folder"); + AliError(Form("Cannot find Run Loader in Default Event Folder")); return; } AliPHOSLoader* gime = dynamic_cast(rn->GetLoader("PHOSLoader")); if (gime == 0x0) { - Error("ExecuteEvent","Can not find PHOS Loader from Run Loader"); + AliError(Form("Cannot find PHOS Loader from Run Loader")); return; } @@ -330,10 +350,11 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) const } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits) +void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits, TVector3 &vInc) { // Calculates the dispersion of the shower at the origine of the RecPoint - + //DP: should we correct dispersion for non-perpendicular hit???????? + Float_t d = 0. ; Float_t wtot = 0. ; @@ -342,7 +363,7 @@ void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = AliPHOSLoader::GetPHOSGeometry(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); // Calculates the center of gravity in the local PHOS-module coordinates @@ -355,13 +376,21 @@ void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits 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 ; + if (fAmp>0 && fEnergyList[iDigit]>0) { + Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; + x += xi * w ; + z += zi * w ; + wtot += w ; + } + else + AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp)); } - x /= wtot ; - z /= wtot ; + if (wtot>0) { + x /= wtot ; + z /= wtot ; + } + else + AliError(Form("Wrong weight %f\n", wtot)); // Calculates the dispersion in coordinates @@ -373,17 +402,26 @@ void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits 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 ; - - + if (fAmp>0 && fEnergyList[iDigit]>0) { + Float_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + d += w*((xi-x)*(xi-x) + (zi-z)*(zi-z) ) ; + wtot+=w ; + } + else + AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp)); } - d /= wtot ; + if (wtot>0) { + d /= wtot ; + } + else + AliError(Form("Wrong weight %f\n", wtot)); + + fDispersion = 0; + if (d>=0) + fDispersion = TMath::Sqrt(d) ; - fDispersion = TMath::Sqrt(d) ; } //______________________________________________________________________________ @@ -393,6 +431,7 @@ void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits // 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% +//DP: non-perpendicular incidence?????????????? Float_t coreRadius = 3 ; @@ -401,7 +440,7 @@ void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = AliPHOSLoader::GetPHOSGeometry(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); Int_t iDigit; @@ -414,13 +453,21 @@ void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits 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 ; + if (fAmp>0 && fEnergyList[iDigit]>0) { + Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; + x += xi * w ; + z += zi * w ; + wtot += w ; + } + else + AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp)); + } + if (wtot>0) { + x /= wtot ; + z /= wtot ; } - x /= wtot ; - z /= wtot ; + else + AliError(Form("Wrong weight %f\n", wtot)); for(iDigit=0; iDigit < fMulDigit; iDigit++) { @@ -435,10 +482,11 @@ void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits fCoreEnergy += fEnergyList[iDigit] ; } + } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) +void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits, TVector3 &vInc) { // Calculates the axis of the shower ellipsoid @@ -451,7 +499,7 @@ void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = AliPHOSLoader::GetPHOSGeometry(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); Int_t iDigit; @@ -463,29 +511,34 @@ void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) Float_t zi ; phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; phosgeom->RelPosInModule(relid, xi, zi); - Double_t w = TMath::Max(0.,logWeight+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 ; + if (fAmp>0 && fEnergyList[iDigit]>0) { + Double_t w = TMath::Max(0.,logWeight+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 ; + } + else + AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp)); } - dxx /= wtot ; - x /= wtot ; - dxx -= x * x ; - dzz /= wtot ; - z /= wtot ; - dzz -= z * z ; - dxz /= wtot ; - dxz -= x * z ; + if (wtot>0) { + dxx /= wtot ; + x /= wtot ; + dxx -= x * x ; + dzz /= wtot ; + z /= wtot ; + dzz -= z * z ; + dxz /= wtot ; + dxz -= x * z ; // //Apply correction due to non-perpendicular incidence // Double_t CosX ; // Double_t CosZ ; // AliPHOSGetter * gime = AliPHOSGetter::Instance() ; // AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); - // Double_t DistanceToIP= (Double_t ) phosgeom->GetIPtoCrystalSurface() ; +// Double_t DistanceToIP= (Double_t ) phosgeom->GetIPtoCrystalSurface() ; // CosX = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+x*x) ; // CosZ = DistanceToIP/TMath::Sqrt(DistanceToIP*DistanceToIP+z*z) ; @@ -495,19 +548,24 @@ void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) // 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[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. ; + 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. ; + } + else { + AliError(Form("Wrong weight %f\n", wtot)); + fLambda[0]=fLambda[1]=0.; + } } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) +void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits, TVector3 &vInc) { // Calculate the shower moments in the eigen reference system // M2x, M2z, M3x, M4z @@ -523,7 +581,7 @@ void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = AliPHOSLoader::GetPHOSGeometry(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; Int_t iDigit; @@ -539,27 +597,38 @@ void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; phosgeom->RelPosInModule(relid, xi, zi); - w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; - x += w * xi ; - z += w * zi ; - dxx += w * xi * xi ; - dzz += w * zi * zi ; - dxz += w * xi * zi ; - wtot += w ; + if (fAmp>0 && fEnergyList[iDigit]>0) { + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + x += w * xi ; + z += w * zi ; + dxx += w * xi * xi ; + dzz += w * zi * zi ; + dxz += w * xi * zi ; + wtot += w ; + } + else + AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp)); + } - x /= wtot ; - z /= wtot ; - dxx /= wtot ; - dzz /= wtot ; - dxz /= wtot ; - dxx -= x * x ; - dzz -= z * z ; - dxz -= x * z ; + if (wtot>0) { + x /= wtot ; + z /= wtot ; + dxx /= wtot ; + dzz /= wtot ; + dxz /= wtot ; + dxx -= x * x ; + dzz -= z * z ; + dxz -= x * z ; // 2) Find covariance matrix eigen values lambda0 and lambda1 - lambda0 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; - lambda1 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; + lambda0 = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; + lambda1 = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; + } + else { + AliError(Form("Wrong weight %f\n", wtot)); + lambda0=lambda1=0.; + } // 3) Find covariance matrix eigen vectors e0 and e1 @@ -596,34 +665,44 @@ void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) phosgeom->RelPosInModule(relid, xiPHOS, ziPHOS); xi = xiPHOS*cosPhi + ziPHOS*sinPhi; zi = ziPHOS*cosPhi - xiPHOS*sinPhi; - w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; - x += w * xi ; - z += w * zi ; - dxx += w * xi * xi ; - dzz += w * zi * zi ; - dxz += w * xi * zi ; - dx3 += w * xi * xi * xi; - dz3 += w * zi * zi * zi ; - dz4 += w * zi * zi * zi * zi ; - wtot += w ; + if (fAmp>0 && fEnergyList[iDigit]>0) { + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + x += w * xi ; + z += w * zi ; + dxx += w * xi * xi ; + dzz += w * zi * zi ; + dxz += w * xi * zi ; + dx3 += w * xi * xi * xi; + dz3 += w * zi * zi * zi ; + dz4 += w * zi * zi * zi * zi ; + wtot += w ; + } + else + AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp)); + } + if (wtot>0) { + x /= wtot ; + z /= wtot ; + dxx /= wtot ; + dzz /= wtot ; + dxz /= wtot ; + dx3 /= wtot ; + dz3 /= wtot ; + dz4 /= wtot ; + dx3 += -3*dxx*x + 2*x*x*x; + dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z; + dxx -= x * x ; + dzz -= z * z ; + dxz -= x * z ; } - x /= wtot ; - z /= wtot ; - dxx /= wtot ; - dzz /= wtot ; - dxz /= wtot ; - dx3 /= wtot ; - dz3 /= wtot ; - dz4 /= wtot ; - dx3 += -3*dxx*x + 2*x*x*x; - dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z; - dxx -= x * x ; - dzz -= z * z ; - dxz -= x * z ; + else + AliError(Form("Wrong weight %f\n", wtot)); // 5) Find an angle between cluster center vector and eigen vector e0 - Float_t phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z)); + Float_t phi = 0; + if ( (x*x+z*z) > 0 ) + phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z)); fM2x = lambda0; fM2z = lambda1; @@ -635,18 +714,22 @@ void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) //____________________________________________________________________________ void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits ) { - // Evaluates all shower parameters - - EvalLocalPosition(logWeight, digits) ; - EvalElipsAxis(logWeight, digits) ; - EvalMoments(logWeight, digits) ; - EvalDispersion(logWeight, digits) ; EvalCoreEnergy(logWeight, digits); EvalTime(digits) ; AliPHOSRecPoint::EvalAll(digits) ; } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits) +void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TVector3 &vtx, TClonesArray * digits ) +{ + // Evaluates all shower parameters + TVector3 vInc ; + EvalLocalPosition(logWeight, vtx, digits,vInc) ; + EvalElipsAxis(logWeight, digits, vInc) ; //they are evaluated with momenta + EvalMoments(logWeight, digits, vInc) ; + EvalDispersion(logWeight, digits, vInc) ; +} +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TVector3 &vtx, TClonesArray * digits, TVector3 &vInc) { // Calculates the center of gravity in the local PHOS-module coordinates Float_t wtot = 0. ; @@ -658,7 +741,7 @@ void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * dig AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = AliPHOSLoader::GetPHOSGeometry(); + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; Int_t iDigit; @@ -669,37 +752,36 @@ void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * dig 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 ; - + if (fAmp>0 && fEnergyList[iDigit]>0) { + Float_t w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp ) ) ; + x += xi * w ; + z += zi * w ; + wtot += w ; + } + else + AliError(Form("Wrong energy %f and/or amplitude %f\n", fEnergyList[iDigit], fAmp)); } - - x /= wtot ; - z /= wtot ; + if (wtot>0) { + x /= wtot ; + z /= wtot ; + } + else + AliError(Form("Wrong weight %f\n", 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]) ; + phosgeom->GetIncidentVector(vtx,GetPHOSMod(),x,z,vInc) ; - //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) ; + Float_t depthx = 0.; + Float_t depthz = 0.; + if (fAmp>0&&vInc.Y()!=0.) { + depthx = ( para * TMath::Log(fAmp) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ; + depthz = ( para * TMath::Log(fAmp) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ; + } + else + AliError(Form("Wrong amplitude %f\n", fAmp)); fLocPos.SetX(x - depthx) ; fLocPos.SetY(0.) ; @@ -763,6 +845,9 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax( AliPHOSDigit ** maxAt, Float_t digit = maxAt[iDigit] ; for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) { + if(iDigit == iDigitN) + continue ; + digitN = (AliPHOSDigit *) digits->At(fDigitsList[iDigitN]) ; if ( AreNeighbours(digit, digitN) ) { @@ -813,8 +898,8 @@ void AliPHOSEmcRecPoint::EvalTime(TClonesArray * digits) 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 * tempo = new Int_t[fMaxDigit]; + Float_t * tempoE = new Float_t[fMaxDigit]; Int_t mult = 0 ; for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){ @@ -828,8 +913,8 @@ void AliPHOSEmcRecPoint::Purify(Float_t threshold){ fMulDigit = mult ; delete [] fDigitsList ; delete [] fEnergyList ; - fDigitsList = new (Int_t[fMulDigit]) ; - fEnergyList = new ( Float_t[fMulDigit]) ; + fDigitsList = new Int_t[fMulDigit]; + fEnergyList = new Float_t[fMulDigit]; fAmp = 0.; //Recalculate total energy for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){ @@ -849,27 +934,27 @@ void AliPHOSEmcRecPoint::Print(Option_t *) const TString message ; message = "AliPHOSEmcRecPoint:\n" ; - message += " digits # = " ; - Info("Print", message.Data()) ; + message += "Digit multiplicity = %d" ; + message += ", cluster Energy = %f" ; + message += ", number of primaries = %d" ; + message += ", rec.point index = %d \n" ; + printf(message.Data(), fMulDigit, fAmp, fMulTrack,GetIndexInList() ) ; Int_t iDigit; + printf(" digits ids = ") ; for(iDigit=0; iDigit