X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSEmcRecPoint.cxx;h=92e5d274bf1ab4523a5b9323551316a4a9642a1c;hb=c2140715d75ba09019bc2373b20c7e2e077062c0;hp=9deeb3d1df3a19548dfabd88c02946d1f4f3f924;hpb=afa51c4ecd8a3a3b5dad7dccfa4cfef0132ace50;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSEmcRecPoint.cxx b/PHOS/AliPHOSEmcRecPoint.cxx index 9deeb3d1df3..92e5d274bf1 100644 --- a/PHOS/AliPHOSEmcRecPoint.cxx +++ b/PHOS/AliPHOSEmcRecPoint.cxx @@ -27,11 +27,10 @@ #include "TH2.h" #include "TMath.h" #include "TCanvas.h" +#include "TGraph.h" // --- Standard library --- -#include - // --- AliRoot header files --- #include "AliGenerator.h" @@ -51,8 +50,10 @@ AliPHOSEmcRecPoint::AliPHOSEmcRecPoint() : AliPHOSRecPoint() fAmp = 0. ; fCoreEnergy = 0 ; fEnergyList = 0 ; + fNExMax = 0 ; //Not unfolded yet fTime = -1. ; fLocPos.SetX(1000000.) ; //Local position should be evaluated + fDebug=0; } @@ -63,13 +64,31 @@ AliPHOSEmcRecPoint::AliPHOSEmcRecPoint(const char * opt) : AliPHOSRecPoint(opt) 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) +{ + // 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() { @@ -480,17 +499,145 @@ void AliPHOSEmcRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) fLambda[1]= 0. ; } +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::EvalMoments(Float_t logWeight,TClonesArray * digits) +{ + // 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 ; + + AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ; + AliPHOSGeometry * phosgeom = (AliPHOSGeometry*)gime->PHOSGeometry(); + + 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); + 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 ; + } + 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 ) ; + + // 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; + 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 ; + } + 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 ; + + // 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)); + + fM2x = lambda0; + fM2z = lambda1; + fM3x = dx3; + fM4z = dz4; + fPhixe = phi; + +} //____________________________________________________________________________ void AliPHOSEmcRecPoint::EvalAll(Float_t logWeight, TClonesArray * digits ) { // Evaluates all shower parameters - AliPHOSRecPoint::EvalAll(logWeight,digits) ; EvalLocalPosition(logWeight, digits) ; EvalElipsAxis(logWeight, digits) ; + EvalMoments(logWeight, digits) ; EvalDispersion(logWeight, digits) ; EvalCoreEnergy(logWeight, digits); EvalTime(digits) ; + AliPHOSRecPoint::EvalAll(logWeight,digits) ; } //____________________________________________________________________________ void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits) @@ -541,16 +688,13 @@ void AliPHOSEmcRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * dig 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) ; - fLocPos.SetX(x - depthx) ; fLocPos.SetY(0.) ; @@ -592,7 +736,7 @@ Int_t AliPHOSEmcRecPoint::GetMultiplicityAtLevel(const Float_t H) const } //____________________________________________________________________________ -Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEnergy, +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 @@ -606,28 +750,28 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEn Int_t iDigit ; for(iDigit = 0; iDigit < fMulDigit; iDigit++) - maxAt[iDigit] = (Int_t) digits->At(fDigitsList[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 *) 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] + locMaxCut) - maxAt[iDigit] = -1 ; + maxAt[iDigit] = 0 ; } else { - maxAt[iDigit] = -1 ; + maxAt[iDigit] = 0 ; // but may be digitN too is not local max ? if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut) - maxAt[iDigitN] = -1 ; + maxAt[iDigitN] = 0 ; } } // if Areneighbours } // while digitN @@ -636,7 +780,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++ ; @@ -645,8 +789,10 @@ Int_t AliPHOSEmcRecPoint::GetNumberOfLocalMax(Int_t * maxAt, Float_t * maxAtEn return iDigitN ; } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::EvalTime(TClonesArray * digits){ - +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++){ @@ -659,32 +805,64 @@ void AliPHOSEmcRecPoint::EvalTime(TClonesArray * digits){ } //____________________________________________________________________________ -void AliPHOSEmcRecPoint::Print(Option_t * option) +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 ; + +} +//____________________________________________________________________________ +void AliPHOSEmcRecPoint::Print(Option_t * option) const { // Print the list of digits belonging to the cluster - cout << "AliPHOSEmcRecPoint: " << endl ; + TString message ; + message = "AliPHOSEmcRecPoint:\n" ; + message += " digits # = " ; + Info("Print", message.Data()) ; Int_t iDigit; - cout << " digits # = " ; for(iDigit=0; iDigit