X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSRecParticle.cxx;h=3bc56ebbf89408a5280eeff4dedf320549396479;hb=aa64ff7a23d3170d300ef54ef2be015bdaa3a3d9;hp=a493f3a799d9351c320796fe2e06d2e2a3210ad4;hpb=2731cd1e25d597f63ec4341cc5b3c4b90ef1d012;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSRecParticle.cxx b/PHOS/AliPHOSRecParticle.cxx index a493f3a799d..3bc56ebbf89 100644 --- a/PHOS/AliPHOSRecParticle.cxx +++ b/PHOS/AliPHOSRecParticle.cxx @@ -27,22 +27,39 @@ // --- Standard library --- -#include // --- AliRoot header files --- - +#include "AliStack.h" +#include "AliPHOSHit.h" +#include "AliPHOSDigit.h" +#include "AliPHOSTrackSegment.h" +#include "AliPHOSEmcRecPoint.h" #include "AliPHOSRecParticle.h" -#include "TPad.h" +#include "AliPHOSLoader.h" +#include "AliPHOSGeometry.h" +#include "AliLog.h" -ClassImp(AliPHOSRecParticle) +//____________________________________________________________________________ +AliPHOSRecParticle::AliPHOSRecParticle(): + fPHOSTrackSegment(0), + fDebug(kFALSE), + fPos() +{ + // ctor + const Int_t nSPECIES = AliPID::kSPECIESN; + for(Int_t i = 0; i(rl->GetLoader("PHOSLoader")); + Int_t emcRPindex = dynamic_cast(phosLoader->TrackSegments()->At(GetPHOSTSIndex()))->GetEmcIndex(); + dynamic_cast(phosLoader->EmcRecPoints()->At(emcRPindex))->GetPrimaries(rv) ; + return rv ; +} + +//____________________________________________________________________________ +const TParticle * AliPHOSRecParticle::GetPrimary() const +{ + // Get the primary particle at the origine of the RecParticle and + // which has deposited the largest energy in SDigits + AliRunLoader* rl = AliRunLoader::GetRunLoader() ; + 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 ; + + // 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 = phosLoader->SDigits() ; + AliPHOSDigit * sdig ; + + for(Int_t i = 0 ; i < sdigits->GetEntriesFast() ; i++){ + sdig = static_cast(sdigits->At(i)) ; + if((sdig->GetId() == AbsId)){ + if((sdig->GetAmp() > amp) && (sdig->GetNprimary())){ + amp = sdig->GetAmp() ; + iPrim = sdig->GetPrimary(1) ; + } + // do not scan rest of list + if(sdig->GetId() > AbsId) + continue ; + } + } + } + if(iPrim >= 0) + return rl->Stack()->Particle(iPrim) ; + else + return 0 ; +} + +//____________________________________________________________________________ +Int_t AliPHOSRecParticle::GetPrimaryIndex() const { - // Retrieves all the primary particles at the origine of this reconstructed particle + // 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 -// AliPHOSTrackSegment * ts = GetPHOSTrackSegment() ; -// Int_t emcnumber = 0 ; -// Int_t * emclist = ts->GetPrimariesEmc(emcnumber) ; + AliRunLoader* rl = AliRunLoader::GetRunLoader() ; + AliPHOSLoader * phosLoader = dynamic_cast(rl->GetLoader("PHOSLoader")); + rl->GetEvent(rl->GetEventNumber()) ; + rl->LoadHits("READ"); + rl->LoadDigits("READ"); + rl->LoadRecPoints("READ"); + rl->LoadTracks("READ"); -// Int_t ppsdlnumber = 0 ; -// Int_t * ppsdllist = ts->GetPrimariesPpsdLow(ppsdlnumber) ; - -// Int_t ppsdunumber = 0 ; -// Int_t * ppsdulist = ts->GetPrimariesPpsdUp(ppsdunumber) ; - -// number = emcnumber + ppsdlnumber + ppsdunumber ; -// Int_t * list = new Int_t[number] ; + // 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; + } -// Int_t index ; -// for ( index = 0 ; index < emcnumber ; index++) -// list[index] = emclist[index] ; - -// Int_t jndex ; -// for ( jndex = 0 ; jndex < ppsdlnumber ; jndex++) { -// assert(index < number) ; -// list[index] = ppsdllist[jndex] ; -// index++ ; -// } - -// for ( jndex = 0 ; jndex < ppsdunumber ; jndex++) { -// assert(index < number) ; -// list[index] = ppsdulist[jndex] ; -// index++ ; -// } - -// delete emclist ; -// delete ppsdllist ; -// delete ppsdulist ; - - return 0 ; //<--- list ; -} + // 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 +{ + // Get one of the primary particles at the origine of the RecParticle + if ( index > GetNPrimariesToRecParticles() ) { + if (fDebug) + Warning("GetPrimary", "AliPHOSRecParticle::GetPrimary -> %d is larger that the number of primaries %d", + index, GetNPrimaries()) ; + return 0 ; + } + else { + Int_t dummy ; + AliRunLoader* rl = AliRunLoader::GetRunLoader() ; + AliPHOSLoader * phosLoader = dynamic_cast(rl->GetLoader("PHOSLoader")); + + 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, Float_t weight) +{ + // Set the probability densities that this reconstructed particle + // has a type of i: + // i particle types + // ---------------------- + // 0 electron + // 1 muon + // 2 pi+- + // 3 K+- + // 4 p/pbar + // 5 photon + // 6 pi0 at high pt + // 7 neutron + // 8 K0L + // 9 Conversion electron + + fPID[type] = weight ; +}