X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSRecParticle.cxx;h=34bc0d47a3de03fa78dac7909d7358a4b39185dc;hb=123109bbc53922afc0faf7d6bf3a62c223f07b3e;hp=0d2d7d5922863fb3937621c6306e0a320a8961a7;hpb=17323043f584137499bf2a6e2f2e2da17448b152;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSRecParticle.cxx b/PHOS/AliPHOSRecParticle.cxx index 0d2d7d59228..34bc0d47a3d 100644 --- a/PHOS/AliPHOSRecParticle.cxx +++ b/PHOS/AliPHOSRecParticle.cxx @@ -29,28 +29,37 @@ // --- AliRoot header files --- +#include "AliStack.h" +#include "AliPHOSHit.h" +#include "AliPHOSDigit.h" +#include "AliPHOSTrackSegment.h" +#include "AliPHOSEmcRecPoint.h" #include "AliPHOSRecParticle.h" -#include "AliPHOSGetter.h" +#include "AliPHOSLoader.h" #include "AliPHOSGeometry.h" +#include "AliLog.h" //____________________________________________________________________________ - AliPHOSRecParticle::AliPHOSRecParticle(): fPHOSTrackSegment(0) , fDebug( kFALSE ) +AliPHOSRecParticle::AliPHOSRecParticle(): + fPHOSTrackSegment(0), + fDebug(kFALSE), + fPos() { // ctor - const Int_t nSPECIES = AliESDtrack::kSPECIESN; + const Int_t nSPECIES = AliPID::kSPECIESN; for(Int_t i = 0; i(gime->TrackSegments()->At(GetPHOSTSIndex()))->GetEmcIndex(); - dynamic_cast(gime->EmcRecPoints()->At(emcRPindex))->GetPrimaries(rv) ; + AliRunLoader* rl = AliRunLoader::Instance() ; + AliPHOSLoader * phosLoader = dynamic_cast(rl->GetLoader("PHOSLoader")); + Int_t emcRPindex = dynamic_cast(phosLoader->TrackSegments()->At(GetPHOSTSIndex()))->GetEmcIndex(); + dynamic_cast(phosLoader->EmcRecPoints()->At(emcRPindex))->GetPrimaries(rv) ; return rv ; } @@ -100,24 +110,32 @@ const TParticle * AliPHOSRecParticle::GetPrimary() const { // Get the primary particle at the origine of the RecParticle and // which has deposited the largest energy in SDigits - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; - if (!gime) - Error("GetPrimary", "Getter not yet instantiated") ; - gime->Event(gime->EventNumber(), "SRTPX") ; + AliRunLoader* rl = AliRunLoader::Instance() ; + AliPHOSLoader * phosLoader = dynamic_cast(rl->GetLoader("PHOSLoader")); + rl->GetEvent(rl->GetEventNumber()) ; + rl->LoadKinematics("READ"); + rl->LoadSDigits("READ"); + if(GetNPrimaries() == 0) return 0 ; if(GetNPrimaries() == 1) return GetPrimary(0) ; Int_t AbsId = 0; Int_t module ; - const AliPHOSGeometry * geom = gime->PHOSGeometry() ; - Double_t x,z ; - geom->ImpactOnEmc(static_cast(Theta()),static_cast(Phi()), module,z,x); + + // Get PHOS Geometry object + AliPHOSGeometry *geom; + if (!(geom = AliPHOSGeometry::GetInstance())) + geom = AliPHOSGeometry::GetInstance("IHEP",""); + Double_t x,z ; +//DP to be fixed: Why do we use this method here??? M.B. use RecPoint??? + Double_t vtx[3]={0.,0.,0.} ; + geom->ImpactOnEmc(vtx,static_cast(Theta()),static_cast(Phi()), module,z,x); Int_t amp = 0 ; Int_t iPrim=-1 ; if(module != 0){ geom->RelPosToAbsId(module,x,z,AbsId) ; - TClonesArray * sdigits = gime->SDigits() ; + TClonesArray * sdigits = phosLoader->SDigits() ; AliPHOSDigit * sdig ; for(Int_t i = 0 ; i < sdigits->GetEntriesFast() ; i++){ @@ -134,11 +152,117 @@ const TParticle * AliPHOSRecParticle::GetPrimary() const } } if(iPrim >= 0) - return gime->Primary(iPrim) ; + return rl->Stack()->Particle(iPrim) ; else return 0 ; } +//____________________________________________________________________________ +Int_t AliPHOSRecParticle::GetPrimaryIndex() const +{ + // Get the primary track index in TreeK which deposits the most energy + // in Digits which forms EmcRecPoint, which produces TrackSegment, + // which the RecParticle is created from + + + AliRunLoader* rl = AliRunLoader::Instance() ; + AliPHOSLoader * phosLoader = dynamic_cast(rl->GetLoader("PHOSLoader")); + rl->GetEvent(rl->GetEventNumber()) ; + rl->LoadHits("READ"); + rl->LoadDigits("READ"); + rl->LoadRecPoints("READ"); + rl->LoadTracks("READ"); + + // Get TrackSegment corresponding to this RecParticle + const AliPHOSTrackSegment *ts = phosLoader->TrackSegment(fPHOSTrackSegment); + + // Get EmcRecPoint corresponding to this TrackSegment + Int_t emcRecPointIndex = ts->GetEmcIndex(); + + const AliPHOSEmcRecPoint *emcRecPoint = phosLoader->EmcRecPoint(emcRecPointIndex); + + // Get the list of digits forming this EmcRecParticle + Int_t nDigits = emcRecPoint->GetDigitsMultiplicity(); + Int_t *digitList = emcRecPoint->GetDigitsList(); + + // Find the digit with maximum amplitude + Int_t maxAmp = 0; + Int_t bestDigitIndex = -1; + for (Int_t iDigit=0; iDigitDigit(digitList[iDigit]); + if (digit->GetAmp() > maxAmp) { + maxAmp = digit->GetAmp(); + bestDigitIndex = iDigit; + } + } + if (bestDigitIndex>-1) { + const AliPHOSDigit * digit = phosLoader->Digit(digitList[bestDigitIndex]); + if (digit==0) return -12345; + } + + // Get the list of primary tracks producing this digit + // and find which track has more track energy. + //Int_t nPrimary = digit->GetNprimary(); + //TParticle *track = 0; + //Double_t energyEM = 0; + //Double_t energyHadron = 0; + //Int_t trackEM = -1; + //Int_t trackHadron = -1; + //for (Int_t iPrim=0; iPrimGetPrimary(iPrim+1); + // if (iPrimary == -1 || TMath::Abs(iPrimary)>10000000) + // continue ; //PH 10000000 is the shift of the track + // //PH indexes in the underlying event + // track = gime->Primary(iPrimary); + // Int_t pdgCode = track->GetPdgCode(); + // Double_t energy = track->Energy(); + // if (pdgCode==22 || TMath::Abs(pdgCode)==11) { + // if (energy > energyEM) { + // energyEM = energy; + // trackEM = iPrimary; + // } + // } + // else { + // if (energy > energyHadron) { + // energyHadron = energy; + // trackHadron = iPrimary; + // } + // } + //} + // Preferences are given to electromagnetic tracks + //if (trackEM != -1) return trackEM; // track is gamma or e+- + //if (trackHadron != -1) return trackHadron; // track is hadron + //return -12345; // no track found :( + + + // Get the list of hits producing this digit, + // find which hit has deposited more energy + // and find the primary track. + + TClonesArray *hits = phosLoader->Hits(); + if (hits==0) return -12345; + + Double_t maxedep = 0; + Int_t maxtrack = -1; + Int_t nHits = hits ->GetEntries(); + Int_t id = (phosLoader->Digit(digitList[bestDigitIndex]))->GetId(); + + for (Int_t iHit=0; iHitHit(iHit); + if(hit->GetId() == id){ + Double_t edep = hit->GetEnergy(); + Int_t track = hit->GetPrimary(); + if(edep > maxedep){ + maxedep = edep; + maxtrack = track; + } + } + } + + if (maxtrack != -1) return maxtrack; // track is hadron + return -12345; // no track found :( +} + //____________________________________________________________________________ const TParticle * AliPHOSRecParticle::GetPrimary(Int_t index) const { @@ -151,17 +275,18 @@ const TParticle * AliPHOSRecParticle::GetPrimary(Int_t index) const } else { Int_t dummy ; - AliPHOSGetter * gime = AliPHOSGetter::Instance() ; + AliRunLoader* rl = AliRunLoader::Instance() ; + AliPHOSLoader * phosLoader = dynamic_cast(rl->GetLoader("PHOSLoader")); - Int_t emcRPindex = dynamic_cast(gime->TrackSegments()->At(GetPHOSTSIndex()))->GetEmcIndex(); - Int_t primaryindex = dynamic_cast(gime->EmcRecPoints()->At(emcRPindex))->GetPrimaries(dummy)[index] ; - return gime->Primary(primaryindex) ; + Int_t emcRPindex = dynamic_cast(phosLoader->TrackSegments()->At(GetPHOSTSIndex()))->GetEmcIndex(); + Int_t primaryindex = dynamic_cast(phosLoader->EmcRecPoints()->At(emcRPindex))->GetPrimaries(dummy)[index] ; + return rl->Stack()->Particle(primaryindex) ; } // return 0 ; } //____________________________________________________________________________ -void AliPHOSRecParticle::SetPID(Int_t type, Double_t weight) +void AliPHOSRecParticle::SetPID(Int_t type, Float_t weight) { // Set the probability densities that this reconstructed particle // has a type of i: