X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSEmcRecPoint.cxx;h=9614579aa4d4e51324cfc40085614a86c5d9dfb8;hb=89557f6dd490ae3be0e3f017b2516d25a16a0099;hp=ac68121bdc5338720df6ed16622cd7835b5ef826;hpb=5830e1d94dc403f63a846006c0f760bfef5db051;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSEmcRecPoint.cxx b/PHOS/AliPHOSEmcRecPoint.cxx index ac68121bdc5..9614579aa4d 100644 --- a/PHOS/AliPHOSEmcRecPoint.cxx +++ b/PHOS/AliPHOSEmcRecPoint.cxx @@ -15,48 +15,88 @@ /* $Id$ */ +/* History of cvs commits: + * + * $Log$ + * 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 -// +//*-- //*-- Author: Dmitri Peressounko (RRC KI & SUBATECH) // --- ROOT system --- -#include "TPad.h" #include "TH2.h" #include "TMath.h" #include "TCanvas.h" +#include "TGraph.h" // --- Standard library --- -#include - // --- AliRoot header files --- - +#include "AliLog.h" +#include "AliPHOSLoader.h" #include "AliGenerator.h" #include "AliPHOSGeometry.h" +#include "AliPHOSDigit.h" #include "AliPHOSEmcRecPoint.h" -#include "AliRun.h" -#include "AliPHOSIndexToObject.h" - + ClassImp(AliPHOSEmcRecPoint) //____________________________________________________________________________ -AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(Float_t W0, Float_t LocMaxCut) - : AliPHOSRecPoint() +AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : + AliPHOSRecPoint(), + fCoreEnergy(0.), fDispersion(0.), + fEnergyList(0), fTime(-1.), fNExMax(0), + fM2x(0.), fM2z(0.), fM3x(0.), fM4z(0.), + fPhixe(0.), fDebug(0) { // ctor + fMulDigit = 0 ; + fAmp = 0. ; + fLocPos.SetX(1000000.) ; //Local position should be evaluated +} +//____________________________________________________________________________ +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.), fDebug(0) +{ + // ctor fMulDigit = 0 ; fAmp = 0. ; - fEnergyList = new Float_t[fMaxDigit]; - fGeom = AliPHOSGeometry::GetInstance() ; - fDelta = ((AliPHOSGeometry *) fGeom)->GetCrystalSize(0) ; - fW0 = W0 ; - fLocMaxCut = LocMaxCut ; fLocPos.SetX(1000000.) ; //Local position should be evaluated - +} + +//____________________________________________________________________________ +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), fDebug(rp.fDebug) +{ + // cpy ctor + fMulDigit = rp.fMulDigit ; + fAmp = rp.fAmp ; + fEnergyList = new Float_t[rp.fMulDigit] ; + Int_t index ; + for(index = 0 ; index < fMulDigit ; index++) + fEnergyList[index] = rp.fEnergyList[index] ; } //____________________________________________________________________________ @@ -71,12 +111,15 @@ AliPHOSEmcRecPoint::~AliPHOSEmcRecPoint() void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy) { // Adds a digit to the RecPoint - // and accumulates the total amplitude and the multiplicity + // 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_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++ ){ @@ -85,10 +128,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] ; @@ -103,16 +146,19 @@ void AliPHOSEmcRecPoint::AddDigit(AliPHOSDigit & digit, Float_t Energy) fEnergyList[fMulDigit] = Energy ; fMulDigit++ ; fAmp += Energy ; + + EvalPHOSMod(&digit) ; } //____________________________________________________________________________ Bool_t AliPHOSEmcRecPoint::AreNeighbours(AliPHOSDigit * digit1, AliPHOSDigit * digit2 ) const { - // Tells if (true) or not (false) two digits are neighbors) + // Tells if (true) or not (false) two digits are neighbors Bool_t aren = kFALSE ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + Int_t relid1[4] ; phosgeom->AbsToRelNumbering(digit1->GetId(), relid1) ; @@ -132,7 +178,10 @@ 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 - + + 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 ; AliPHOSEmcRecPoint * clu = (AliPHOSEmcRecPoint *)obj ; @@ -142,16 +191,16 @@ Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const 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 @@ -167,20 +216,18 @@ Int_t AliPHOSEmcRecPoint::Compare(const TObject * obj) const return rv ; } - //______________________________________________________________________________ -void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) +void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t, Int_t) /*const*/ { + // 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; - - AliPHOSIndexToObject * please = AliPHOSIndexToObject::GetInstance() ; + + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); static TGraph * digitgraph = 0 ; @@ -188,12 +235,29 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) TH2F * histo = 0 ; TCanvas * histocanvas ; + + + //try to get run loader from default event folder + AliRunLoader* rn = AliRunLoader::GetRunLoader(AliConfig::GetDefaultEventFolderName()); + if (rn == 0x0) + { + AliError(Form("Cannot find Run Loader in Default Event Folder")); + return; + } + AliPHOSLoader* gime = dynamic_cast(rn->GetLoader("PHOSLoader")); + if (gime == 0x0) + { + AliError(Form("Cannot find PHOS Loader from Run Loader")); + return; + } + + + const TClonesArray * digits = gime->Digits() ; switch (event) { case kButton1Down: { AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; Int_t iDigit; Int_t relid[4] ; @@ -209,7 +273,7 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) Float_t zimin = 999. ; for(iDigit=0; iDigitGimeDigit(fDigitsList[iDigit]) ) ; + digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; phosgeom->RelPosInModule(relid, xi[iDigit], zi[iDigit]); if ( xi[iDigit] > ximax ) @@ -241,7 +305,7 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) Float_t x, z ; for(iDigit=0; iDigitGimeDigit(fDigitsList[iDigit]) ) ; + digit = (AliPHOSDigit *) digits->At(fDigitsList[iDigit]) ; phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; phosgeom->RelPosInModule(relid, x, z); histo->Fill(x, z, fEnergyList[iDigit] ) ; @@ -255,8 +319,8 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) digitgraph-> Paint("P") ; } - Print() ; - histocanvas = new TCanvas("cluser", "a single cluster", 600, 500) ; + // Print() ; + histocanvas = new TCanvas("cluster", "a single cluster", 600, 500) ; histocanvas->Draw() ; histo->Draw("lego1") ; @@ -277,64 +341,125 @@ void AliPHOSEmcRecPoint::ExecuteEvent(Int_t event, Int_t px, Int_t py) } //____________________________________________________________________________ -Float_t AliPHOSEmcRecPoint::GetDispersion() const +void AliPHOSEmcRecPoint::EvalDispersion(Float_t logWeight,TClonesArray * digits) { // Calculates the dispersion of the shower at the origine of the RecPoint - AliPHOSIndexToObject * please = AliPHOSIndexToObject::GetInstance() ; - - Float_t d = 0 ; - Float_t wtot = 0 ; + Float_t d = 0. ; + Float_t wtot = 0. ; - TVector3 locpos; - GetLocalPosition(locpos); - Float_t x = locpos.X() ; - Float_t z = locpos.Z() ; + Float_t x = 0.; + Float_t z = 0.; AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); + + // Calculates the center of gravity in the local PHOS-module coordinates 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); + 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 ; + } + else + AliError(Form("Wrong weight %f\n", wtot)); + + +// Calculates the dispersion in coordinates + wtot = 0.; for(iDigit=0; iDigit < fMulDigit; iDigit++) { - digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[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.,fW0+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)); - return TMath::Sqrt(d) ; + fDispersion = 0; + if (d>=0) + fDispersion = TMath::Sqrt(d) ; + } //______________________________________________________________________________ -Float_t AliPHOSEmcRecPoint::CoreEnergy() +void AliPHOSEmcRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) { - //This function calculates energy in the core, - //i.e. within radius rad = 3cm. Beyond this radius - //in accoradnce with shower profile energy deposition + // 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% - AliPHOSIndexToObject * please = AliPHOSIndexToObject::GetInstance() ; - - Float_t eCore = 0 ; Float_t coreRadius = 3 ; - TVector3 locpos; - GetLocalPosition(locpos); - Float_t x = locpos.X() ; - Float_t z = locpos.Z() ; + Float_t x = 0 ; + Float_t z = 0 ; AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; - + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); + 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); + 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 ; + } + else + AliError(Form("Wrong weight %f\n", wtot)); + + for(iDigit=0; iDigit < fMulDigit; iDigit++) { - digit = (AliPHOSDigit *) ( please->GimeDigit(fDigitsList[iDigit]) ) ; + digit = (AliPHOSDigit *) ( digits->At(fDigitsList[iDigit]) ) ; Int_t relid[4] ; Float_t xi ; Float_t zi ; @@ -342,19 +467,16 @@ Float_t AliPHOSEmcRecPoint::CoreEnergy() phosgeom->RelPosInModule(relid, xi, zi); Float_t distance = TMath::Sqrt((xi-x)*(xi-x)+(zi-z)*(zi-z)) ; if(distance < coreRadius) - eCore += fEnergyList[iDigit] ; + fCoreEnergy += fEnergyList[iDigit] ; } -return eCore ; } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda) +void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) { // Calculates the axis of the shower ellipsoid - AliPHOSIndexToObject * please = AliPHOSIndexToObject::GetInstance() ; - Double_t wtot = 0. ; Double_t x = 0.; Double_t z = 0.; @@ -363,37 +485,47 @@ void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda) Double_t dxz = 0.; AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance(); + Int_t iDigit; + for(iDigit=0; iDigitGimeDigit(fDigitsList[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); - Double_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 ; + 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 ; -// Double_t DistanceToIP= (Double_t ) ((AliPHOSGeometry *) fGeom)->GetIPtoCrystalSurface() ; +// AliPHOSGetter * gime = AliPHOSGetter::Instance() ; +// 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) ; @@ -403,30 +535,185 @@ void AliPHOSEmcRecPoint::GetElipsAxis(Float_t * lambda) // dxz = dxz/(CosX*CosZ) ; - lambda[0] = 0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; - if(lambda[0] > 0) - lambda[0] = TMath::Sqrt(lambda[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]) ; - lambda[1] = 0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ) ; - if(lambda[1] > 0) //To avoid exception if numerical errors lead to negative lambda. - lambda[1] = TMath::Sqrt(lambda[1]) ; - else - lambda[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::EvalAll(void ) +void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) { - AliPHOSRecPoint::EvalAll() ; - EvalLocalPosition() ; + // Calculate the shower moments in the eigen reference system + // M2x, M2z, M3x, M4z + // Calculate the angle between the shower position vector and the eigen vector + + Double_t wtot = 0. ; + Double_t x = 0.; + Double_t z = 0.; + Double_t dxx = 0.; + Double_t dzz = 0.; + Double_t dxz = 0.; + Double_t lambda0=0, lambda1=0; + + AliPHOSDigit * digit ; + + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; + + Int_t iDigit; + + // 1) Find covariance matrix elements: + // || dxx dxz || + // || dxz dzz || + + Float_t xi ; + Float_t zi ; + Int_t relid[4] ; + Double_t w; + for(iDigit=0; iDigitAt(fDigitsList[iDigit]) ; + phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; + phosgeom->RelPosInModule(relid, xi, zi); + 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)); + + } + 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 ) ; + } + else { + AliError(Form("Wrong weight %f\n", wtot)); + lambda0=lambda1=0.; + } + + // 3) Find covariance matrix eigen vectors e0 and e1 + + TVector2 e0,e1; + if (dxz != 0) + e0.Set(1.,(lambda0-dxx)/dxz); + else + e0.Set(0.,1.); + + e0 = e0.Unit(); + e1.Set(-e0.Y(),e0.X()); + + // 4) Rotate cluster tensor from (x,z) to (e0,e1) system + // and calculate moments M3x and M4z + + Float_t cosPhi = e0.X(); + Float_t sinPhi = e0.Y(); + + Float_t xiPHOS ; + Float_t ziPHOS ; + Double_t dx3, dz3, dz4; + wtot = 0.; + x = 0.; + z = 0.; + dxx = 0.; + dzz = 0.; + dxz = 0.; + dx3 = 0.; + dz3 = 0.; + dz4 = 0.; + for(iDigit=0; iDigitAt(fDigitsList[iDigit]) ; + phosgeom->AbsToRelNumbering(digit->GetId(), relid) ; + phosgeom->RelPosInModule(relid, xiPHOS, ziPHOS); + xi = xiPHOS*cosPhi + ziPHOS*sinPhi; + zi = ziPHOS*cosPhi - xiPHOS*sinPhi; + 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 ; + } + else + AliError(Form("Wrong weight %f\n", wtot)); + + // 5) Find an angle between cluster center vector and eigen vector e0 + + 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; + fM3x = dx3; + fM4z = dz4; + fPhixe = phi; + } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::EvalLocalPosition(void ) +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) { // Calculates the center of gravity in the local PHOS-module coordinates - - AliPHOSIndexToObject * please = AliPHOSIndexToObject::GetInstance() ; - Float_t wtot = 0. ; Int_t relid[4] ; @@ -436,34 +723,44 @@ void AliPHOSEmcRecPoint::EvalLocalPosition(void ) AliPHOSDigit * digit ; - AliPHOSGeometry * phosgeom = (AliPHOSGeometry *) fGeom ; + AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ; Int_t iDigit; for(iDigit=0; iDigitGimeDigit(fDigitsList[iDigit]) ); + digit = (AliPHOSDigit *) digits->At(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 ; - + 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) ; - + //We should check all 3 possibilities to avoid seg.v. + if(gAlice && gAlice->GetMCApp() && gAlice->Generator()) + gAlice->Generator()->GetOrigin(xo,yo,zo) ; + else{ + xo=yo=zo=0.; + } Float_t phi = phosgeom->GetPHOSAngle(relid[0]) ; //Transform to the local ref.frame @@ -471,21 +768,25 @@ void AliPHOSEmcRecPoint::EvalLocalPosition(void ) xoL = xo*TMath::Cos(phi)-yo*TMath::Sin(phi) ; yoL = xo*TMath::Sin(phi)+yo*TMath::Cos(phi) ; - Float_t radius = TMath::Sqrt((xoL-x)*(xoL-x)+ - (phosgeom->GetIPtoCrystalSurface()-yoL)*(phosgeom->GetIPtoCrystalSurface()-yoL)+ - (zo-z)*(zo-z)); + 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) { + depthx = ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencephi) ; + depthz = ( para * TMath::Log(fAmp) + parb ) * TMath::Sin(incidencetheta) ; + } + else + AliError(Form("Wrong amplitude %f\n", fAmp)); fLocPos.SetX(x - depthx) ; fLocPos.SetY(0.) ; fLocPos.SetZ(z - depthz) ; + fLocPosM = 0 ; } //____________________________________________________________________________ @@ -506,7 +807,7 @@ Float_t AliPHOSEmcRecPoint::GetMaximalEnergy(void) const } //____________________________________________________________________________ -Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) const +Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(Float_t H) const { // Calculates the multiplicity of digits with energy larger than H*energy @@ -521,12 +822,11 @@ Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) const } //____________________________________________________________________________ -Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEnergy) const +Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax( AliPHOSDigit ** maxAt, Float_t * maxAtEnergy, + Float_t locMaxCut,TClonesArray * digits) const { // Calculates the number of local maxima in the cluster using fLocalMaxCut as the minimum - // energy difference between two local maxima - - AliPHOSIndexToObject * please = AliPHOSIndexToObject::GetInstance() ; + // energy difference between two local maxima AliPHOSDigit * digit ; AliPHOSDigit * digitN ; @@ -535,29 +835,32 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEn Int_t iDigitN ; Int_t iDigit ; - for(iDigit = 0; iDigit < fMulDigit; iDigit++){ - maxAt[iDigit] = (Int_t) ( please->GimeDigit(fDigitsList[iDigit]) ) ; - } + for(iDigit = 0; iDigit < fMulDigit; iDigit++) + maxAt[iDigit] = (AliPHOSDigit*) digits->At(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 *) ( please->GimeDigit(fDigitsList[iDigitN]) ) ; + if(iDigit == iDigitN) + continue ; + + 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 @@ -566,7 +869,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++ ; @@ -574,94 +877,83 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEn } return iDigitN ; } +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::EvalTime(TClonesArray * digits) +{ + // Define a rec.point time as a time in the cell with the maximum energy + + 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; + } + } + fTime = ((AliPHOSDigit*) digits->At(fDigitsList[maxAt]))->GetTime() ; + +} +//____________________________________________________________________________ +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]; + + fAmp = 0.; //Recalculate total energy + for(Int_t iDigit=0;iDigit< fMulDigit ;iDigit++){ + fDigitsList[iDigit] = tempo[iDigit]; + fEnergyList[iDigit] = tempoE[iDigit] ; + fAmp+=tempoE[iDigit]; + } + + delete [] tempo ; + delete [] tempoE ; - -// //____________________________________________________________________________ -// AliPHOSEmcRecPoint& AliPHOSEmcRecPoint::operator = (AliPHOSEmcRecPoint Clu) -// { -// int * dl = Clu.GetDigitsList() ; - -// if(fDigitsList) -// delete fDigitsList ; - -// AliPHOSDigit * digit ; - -// Int_t iDigit; - -// for(iDigit=0; iDigitGimeDigit( 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 ; - cout << " Stored at position " << GetIndexInList() << 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; -// fEnergyList = new Float_t[fMulDigit] ; -// R__b.ReadFastArray(fEnergyList, fMulDigit); -// R__b >> fLocMaxCut; -// R__b >> fW0; -// } else { -// R__b.WriteVersion(AliPHOSEmcRecPoint::IsA()); -// AliPHOSRecPoint::Streamer(R__b); -// R__b << fDelta; -// R__b.WriteFastArray(fEnergyList, fMulDigit); -// R__b << fLocMaxCut; -// R__b << fW0; -// } -// } +