X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PWGJE%2FEMCALJetTasks%2FAliJetModelBaseTask.cxx;h=abdb744e87439c90af8581d6d3f7d55cdfc5f5df;hb=eb556d0c8c2c423e23497d5ddf7fb2e094decefc;hp=13025ee6a88feb5fa2b9a76abe9c68d56981e4d3;hpb=60d77596e010f48971b161d70226cc23666070df;p=u%2Fmrichter%2FAliRoot.git diff --git a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx index 13025ee6a88..abdb744e874 100644 --- a/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx +++ b/PWGJE/EMCALJetTasks/AliJetModelBaseTask.cxx @@ -11,6 +11,7 @@ #include #include #include +#include #include "AliVEvent.h" #include "AliAODCaloCluster.h" @@ -25,6 +26,7 @@ #include "AliPicoTrack.h" #include "AliEMCALGeometry.h" #include "AliLog.h" +#include "AliNamedArrayI.h" ClassImp(AliJetModelBaseTask) @@ -40,6 +42,7 @@ AliJetModelBaseTask::AliJetModelBaseTask() : fOutCellsName(), fMCParticlesName(), fOutMCParticlesName(), + fIsMC(kFALSE), fSuffix(), fEtaMin(-1), fEtaMax(1), @@ -53,7 +56,13 @@ AliJetModelBaseTask::AliJetModelBaseTask() : fNTracks(0), fMarkMC(99999), fPtSpectrum(0), + fPtPhiEvPlDistribution(0), + fDensitySpectrum(0), + fDifferentialV2(0), + fAddV2(kFALSE), + fFlowFluctuations(kFALSE), fQAhistos(kFALSE), + fPsi(0), fIsInit(0), fGeom(0), fClusters(0), @@ -72,6 +81,10 @@ AliJetModelBaseTask::AliJetModelBaseTask() : fOutput(0) { // Default constructor. + + fVertex[0] = 0; + fVertex[1] = 0; + fVertex[2] = 0; } //________________________________________________________________________ @@ -86,6 +99,7 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) : fOutCellsName(""), fMCParticlesName(""), fOutMCParticlesName(""), + fIsMC(kFALSE), fSuffix("Processed"), fEtaMin(-1), fEtaMax(1), @@ -99,7 +113,13 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) : fNTracks(1), fMarkMC(99999), fPtSpectrum(0), + fPtPhiEvPlDistribution(0), + fDensitySpectrum(0), + fDifferentialV2(0), + fAddV2(kFALSE), + fFlowFluctuations(kFALSE), fQAhistos(drawqa), + fPsi(0), fIsInit(0), fGeom(0), fClusters(0), @@ -122,6 +142,10 @@ AliJetModelBaseTask::AliJetModelBaseTask(const char *name, Bool_t drawqa) : if (fQAhistos) { DefineOutput(1, TList::Class()); } + + fVertex[0] = 0; + fVertex[1] = 0; + fVertex[2] = 0; } //________________________________________________________________________ @@ -155,6 +179,14 @@ void AliJetModelBaseTask::UserExec(Option_t *) if (!fIsInit) return; + fVertex[0] = 0; + fVertex[1] = 0; + fVertex[2] = 0; + + const AliVVertex *vert = InputEvent()->GetPrimaryVertex(); + if (vert) + vert->GetXYZ(fVertex); + if (fCopyArray) { if (fOutTracks) fOutTracks->Delete(); @@ -164,6 +196,16 @@ void AliJetModelBaseTask::UserExec(Option_t *) fOutMCParticles->Delete(); } + if (fDensitySpectrum) { + fNTracks = TMath::Nint(fDensitySpectrum->GetRandom()); + fNCells = TMath::Nint(fDensitySpectrum->GetRandom()); + fNClusters = TMath::Nint(fDensitySpectrum->GetRandom()); + } + + // Clear map + if (fOutMCParticlesMap) + fOutMCParticlesMap->Clear(); + AliVCaloCells *tempCaloCells = 0; if (fCaloCells) { @@ -174,9 +216,12 @@ void AliJetModelBaseTask::UserExec(Option_t *) } } + if (fPtPhiEvPlDistribution || fAddV2) + fPsi = gRandom->Rndm() * TMath::Pi(); + Run(); - if (fCaloCells) { + if (fCaloCells && !fCopyArray) { delete fCaloCells; fCaloCells = tempCaloCells; } @@ -242,7 +287,7 @@ Bool_t AliJetModelBaseTask::ExecOnce() } } - if (fNTracks > 0 && !fTracksName.IsNull()) { + if (!fTracksName.IsNull()) { fTracks = dynamic_cast(InputEvent()->FindListObject(fTracksName)); if (!fTracks) { AliWarning(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data())); @@ -274,7 +319,11 @@ Bool_t AliJetModelBaseTask::ExecOnce() } } - if (fNClusters > 0 && !fCaloName.IsNull()) { + if(fAddV2 && (!fDifferentialV2)) { + AliWarning(Form("%s: Cannot add v2 without diffential v2!", GetName())); + } + + if (!fCaloName.IsNull()) { fClusters = dynamic_cast(InputEvent()->FindListObject(fCaloName)); if (!fClusters) { @@ -290,8 +339,16 @@ Bool_t AliJetModelBaseTask::ExecOnce() fOutCaloName = fCaloName; if (fCopyArray) fOutCaloName += fSuffix; + TString className; + if (fClusters) + className = fClusters->GetClass()->GetName(); + else if (fEsdMode) + className = "AliESDCaloCluster"; + else + className = "AliAODCaloCluster"; + if (fCopyArray || !fClusters) { - fOutClusters = new TClonesArray(fClusters->GetClass()->GetName(), fClusters->GetSize()); + fOutClusters = new TClonesArray(className.Data()); fOutClusters->SetName(fOutCaloName); if (InputEvent()->FindListObject(fOutCaloName)) { AliFatal(Form("%s: Collection %s is already present in the event!", GetName(), fOutCaloName.Data())); @@ -319,13 +376,13 @@ Bool_t AliJetModelBaseTask::ExecOnce() return kFALSE; } - fMCParticlesMap = dynamic_cast(InputEvent()->FindListObject(fMCParticlesName + "_Map")); + fMCParticlesMap = dynamic_cast(InputEvent()->FindListObject(fMCParticlesName + "_Map")); if (!fMCParticlesMap) { AliWarning(Form("%s: Could not retrieve map for MC particles %s! Will assume MC labels consistent with indexes...", GetName(), fMCParticlesName.Data())); - fMCParticlesMap = new TH1I(fMCParticlesName + "_Map", fMCParticlesName + "_Map",9999,0,1); - for (Int_t i = 0; i < 9999; i++) { - fMCParticlesMap->SetBinContent(i,i); + fMCParticlesMap = new AliNamedArrayI(fMCParticlesName + "_Map", 99999); + for (Int_t i = 0; i < 99999; i++) { + fMCParticlesMap->AddAt(i,i); } } } @@ -346,7 +403,7 @@ Bool_t AliJetModelBaseTask::ExecOnce() InputEvent()->AddObject(fOutMCParticles); } - fOutMCParticlesMap = new TH1I(fOutMCParticlesName + "_Map", fOutMCParticlesName + "_Map",9999,0,1); + fOutMCParticlesMap = new AliNamedArrayI(fOutMCParticlesName + "_Map",99999); if (InputEvent()->FindListObject(fOutMCParticlesName + "_Map")) { AliFatal(Form("%s: Map %s_Map is already present in the event!", GetName(), fOutMCParticlesName.Data())); return kFALSE; @@ -362,37 +419,20 @@ Bool_t AliJetModelBaseTask::ExecOnce() } } - if (!fCaloName.IsNull() || !fCellsName.IsNull()) { - if (!fGeom) { - if (fGeomName.Length() > 0) { - fGeom = AliEMCALGeometry::GetInstance(fGeomName); - if (!fGeom) { - AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data())); - return kFALSE; - } - } else { - fGeom = AliEMCALGeometry::GetInstance(); - if (!fGeom) { - AliFatal("Could not get default geometry!"); - return kFALSE; - } + if (!fGeom) { + if (fGeomName.Length() > 0) { + fGeom = AliEMCALGeometry::GetInstance(fGeomName); + if (!fGeom) { + AliFatal(Form("Could not get geometry with name %s!", fGeomName.Data())); + return kFALSE; + } + } else { + fGeom = AliEMCALGeometry::GetInstance(); + if (!fGeom) { + AliFatal("Could not get default geometry!"); + return kFALSE; } } - - const Double_t EmcalMinEta = fGeom->GetArm1EtaMin(); - const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax(); - const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad(); - const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad(); - - if (fEtaMax > EmcalMaxEta) fEtaMax = EmcalMaxEta; - if (fEtaMax < EmcalMinEta) fEtaMax = EmcalMinEta; - if (fEtaMin > EmcalMaxEta) fEtaMin = EmcalMaxEta; - if (fEtaMin < EmcalMinEta) fEtaMin = EmcalMinEta; - - if (fPhiMax > EmcalMaxPhi) fPhiMax = EmcalMaxPhi; - if (fPhiMax < EmcalMinPhi) fPhiMax = EmcalMinPhi; - if (fPhiMin > EmcalMaxPhi) fPhiMin = EmcalMaxPhi; - if (fPhiMin < EmcalMinPhi) fPhiMin = EmcalMinPhi; } return kTRUE; @@ -414,28 +454,6 @@ Int_t AliJetModelBaseTask::SetNumberOfOutCells(Int_t n) return n; } -//________________________________________________________________________ -void AliJetModelBaseTask::CopyCells() -{ - if (!fCaloCells) - return; - - for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) { - Int_t mclabel = 0; - Double_t efrac = 0.; - Double_t time = -1; - Short_t cellNum = -1; - Double_t amp = -1; - - fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac); - fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac); - } - - fAddedCells = fCaloCells->GetNumberOfCells(); - - AliDebug(2, Form("%d cells from the current event", fAddedCells)); -} - //________________________________________________________________________ Int_t AliJetModelBaseTask::AddCell(Double_t e, Double_t eta, Double_t phi) { @@ -471,22 +489,107 @@ Int_t AliJetModelBaseTask::AddCell(Double_t e, Int_t absId, Double_t time, Int_t { // Add a cell to the event. - if (label == 0) - label = fMarkMC; - else - label += fMCLabelShift; + if (label < 0) + label = 0; - Bool_t r = fOutCaloCells->SetCell(fAddedCells, absId, e, time, label, 0); + label += fMarkMC + fMCLabelShift; + + Short_t pos = -1; + if (fCaloCells) + pos = fCaloCells->GetCellPosition(absId); + + Double_t efrac = 1; + Bool_t increaseOnSuccess = kFALSE; + + if (pos < 0) { + increaseOnSuccess = kTRUE; + pos = fAddedCells; + } + else { + Short_t cellNumber = -1; + Double_t old_e = 0; + Double_t old_time = 0; + Int_t old_label = 0; + Double_t old_efrac = 0; + fOutCaloCells->GetCell(pos, cellNumber, old_e, old_time, old_label, old_efrac); + + efrac = e / (old_e + e); + + if (old_label > 0 && e < old_efrac * old_e) { + label = old_label; + efrac = old_efrac; + time = old_time; + } + + e += old_e; + } + + Bool_t r = fOutCaloCells->SetCell(pos, absId, e, time, label, efrac); if (r) { - fAddedCells++; + if (increaseOnSuccess) + fAddedCells++; return fAddedCells; } - else { + else return 0; - } } +//________________________________________________________________________ +AliVCluster* AliJetModelBaseTask::AddCluster(AliVCluster *oc) +{ + // Add a cluster to the event. + + const Int_t nClusters = fOutClusters->GetEntriesFast(); + AliVCluster *dc = static_cast(fOutClusters->New(nClusters)); + dc->SetType(AliVCluster::kEMCALClusterv1); + dc->SetE(oc->E()); + Float_t pos[3] = {0}; + oc->GetPosition(pos); + dc->SetPosition(pos); + dc->SetNCells(oc->GetNCells()); + dc->SetCellsAbsId(oc->GetCellsAbsId()); + dc->SetCellsAmplitudeFraction(oc->GetCellsAmplitudeFraction()); + dc->SetID(oc->GetID()); + dc->SetDispersion(oc->GetDispersion()); + dc->SetEmcCpvDistance(-1); + dc->SetChi2(-1); + dc->SetTOF(oc->GetTOF()); //time-of-flight + dc->SetNExMax(oc->GetNExMax()); //number of local maxima + dc->SetM02(oc->GetM02()); + dc->SetM20(oc->GetM20()); + dc->SetDistanceToBadChannel(oc->GetDistanceToBadChannel()); + + //MC + UInt_t nlabels = oc->GetNLabels(); + Int_t *labels = oc->GetLabels(); + TArrayI parents; + if (nlabels != 0 && labels) + parents.Set(nlabels, labels); + else { + nlabels = 1; + parents.Set(1); + parents[0] = 0; + } + + if (fMarkMC+fMCLabelShift != 0) { + for (UInt_t i = 0; i < nlabels; i++) { + parents[i] += fMarkMC+fMCLabelShift; + } + } + + AliESDCaloCluster *esdClus = dynamic_cast(dc); + if (esdClus) { + esdClus->AddLabels(parents); + } + else { + AliAODCaloCluster *aodClus = dynamic_cast(dc); + if (aodClus) + aodClus->SetLabel(parents.GetArray(), nlabels); + } + + return dc; +} //________________________________________________________________________ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Double_t eta, Double_t phi, Int_t label) @@ -555,10 +658,10 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t labe cluster->SetEmcCpvDistance(-1); //MC label - if (label == 0) - label = fMarkMC; - else - label += fMCLabelShift; + if (label < 0) + label = 0; + + label += fMarkMC+fMCLabelShift; if (fEsdMode) { AliESDCaloCluster *esdClus = static_cast(cluster); @@ -574,33 +677,42 @@ AliVCluster* AliJetModelBaseTask::AddCluster(Double_t e, Int_t absId, Int_t labe } //________________________________________________________________________ -AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Bool_t ise, Int_t label) +AliPicoTrack* AliJetModelBaseTask::AddTrack(Double_t pt, Double_t eta, Double_t phi, Byte_t type, Double_t etaemc, Double_t phiemc, Double_t ptemc, Bool_t ise, Int_t label, Short_t charge, Double_t mass) { // Add a track to the event. - - const Int_t nTracks = fOutTracks->GetEntriesFast(); - if (pt < 0) - pt = GetRandomPt(); - if (eta < -100) - eta = GetRandomEta(); - if (phi < 0) - phi = GetRandomPhi(); + if (pt < 0 && eta < -100 && phi < 0) { + GetRandomParticle(pt,eta,phi); + } + else { + if (pt < 0) + pt = GetRandomPt(); + if (eta < -100) + eta = GetRandomEta(); + if (phi < 0) + phi = GetRandomPhi(pt); + } + + if (label >= 0) + label += fMarkMC+fMCLabelShift; + else if (label < 0) + label -= fMarkMC+fMCLabelShift; + + if(fAddV2) AddV2(phi, pt); - if (label == 0) - label = fMarkMC; - else - label += fMCLabelShift; + const Int_t nTracks = fOutTracks->GetEntriesFast(); AliPicoTrack *track = new ((*fOutTracks)[nTracks]) AliPicoTrack(pt, eta, phi, - 1, + charge, label, type, etaemc, - phiemc, - ise); + phiemc, + ptemc, + ise, + mass); return track; } @@ -612,11 +724,61 @@ AliAODMCParticle* AliJetModelBaseTask::AddMCParticle(AliAODMCParticle *part, Int AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part); - fOutMCParticlesMap->SetBinContent(origIndex + fMCLabelShift, nPart); + if (origIndex + fMCLabelShift >= fOutMCParticlesMap->GetSize()) + fOutMCParticlesMap->Set((origIndex + fMCLabelShift)*2); + + fOutMCParticlesMap->AddAt(nPart, origIndex + fMCLabelShift); + AliDebug(2, Form("Setting bin %d to %d (fMCLabelShift=%d, origIndex=%d)", + origIndex + fMCLabelShift, fOutMCParticlesMap->At(origIndex + fMCLabelShift), fMCLabelShift, origIndex)); return aodpart; } +//_____________________________________________________________________________ +void AliJetModelBaseTask::AddV2(Double_t &phi, Double_t &pt) const +{ + // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations + Double_t phi0(phi), v2(0.), f(0.), fp(0.), phiprev(0.); + if(fDifferentialV2) v2 = fDifferentialV2->Eval(pt); + if(TMath::AreEqualAbs(v2, 0, 1e-5)) return; + // introduce flow fluctuations (gaussian) + if(fFlowFluctuations) v2 += TMath::Sqrt(2*(v2*.25)*(v2*.25))*TMath::ErfInverse(2*(gRandom->Uniform(0, fFlowFluctuations))-1); + for (Int_t i(0); i < 100; i++) { + phiprev=phi; //store last value for comparison + f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi)); + fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi)); //first derivative + phi -= f/fp; + if (TMath::AreEqualAbs(phiprev, phi, 1e-10)) break; + } +} + +//_____________________________________________________________________________ +void AliJetModelBaseTask::CopyCells() +{ + if (!fCaloCells) + return; + + fAddedCells = 0; + fCaloCells->Sort(); + for (Short_t i = 0; i < fCaloCells->GetNumberOfCells(); i++) { + Int_t mclabel = 0; + Double_t efrac = 0.; + Double_t time = -1; + Short_t cellNum = -1; + Double_t amp = -1; + + fCaloCells->GetCell(i, cellNum, amp, time, mclabel, efrac); + + if (!fIsMC) + mclabel = 0; + + fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac); + fAddedCells++; + } + + AliDebug(2, Form("%d cells from the current event", fAddedCells)); +} + //________________________________________________________________________ void AliJetModelBaseTask::CopyClusters() { @@ -632,7 +794,12 @@ void AliJetModelBaseTask::CopyClusters() AliESDCaloCluster *esdcluster = static_cast(fClusters->At(i)); if (!esdcluster || !esdcluster->IsEMCAL()) continue; - new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster); + AliESDCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliESDCaloCluster(*esdcluster); + if (!fIsMC) { + TArrayI *labels = clus->GetLabelsArray(); + if (labels) + labels->Reset(); + } nCopiedClusters++; } } @@ -641,7 +808,9 @@ void AliJetModelBaseTask::CopyClusters() AliAODCaloCluster *aodcluster = static_cast(fClusters->At(i)); if (!aodcluster || !aodcluster->IsEMCAL()) continue; - new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster); + AliAODCaloCluster *clus = new ((*fOutClusters)[nCopiedClusters]) AliAODCaloCluster(*aodcluster); + if (!fIsMC) + clus->SetLabel(0,0); nCopiedClusters++; } } @@ -658,10 +827,12 @@ void AliJetModelBaseTask::CopyTracks() const Int_t nTracks = fTracks->GetEntriesFast(); Int_t nCopiedTracks = 0; for (Int_t i = 0; i < nTracks; ++i) { - AliPicoTrack *track = static_cast(fTracks->At(i)); - if (!track) + AliPicoTrack *picotrack = static_cast(fTracks->At(i)); + if (!picotrack) continue; - new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*track); + AliPicoTrack *track = new ((*fOutTracks)[nCopiedTracks]) AliPicoTrack(*picotrack); + if (!fIsMC && track->GetLabel() != 0) + track->SetLabel(0); nCopiedTracks++; } } @@ -685,18 +856,19 @@ void AliJetModelBaseTask::CopyMCParticles() nCopiedPart++; } - if (!fMCParticlesMap) + if (!fMCParticlesMap || !fOutMCParticlesMap) return; - Int_t shift = 0; + if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize()) + fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2); - for (Int_t i = 0; i < fMCParticlesMap->GetNbinsX()+2; i++) { - fOutMCParticlesMap->SetBinContent(i, fMCParticlesMap->GetBinContent(i)); - if (fMCParticlesMap->GetBinContent(i) != 0) - shift = i; + for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) { + fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i); + if (fMCParticlesMap->At(i) >= 0) + fMCLabelShift = i; } - fMCLabelShift = shift; + AliDebug(2,Form("MC particles copied. fMCLabelShift=%d",fMCLabelShift)); } //________________________________________________________________________ @@ -709,9 +881,9 @@ void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &abs Double_t rndPhi = phi; do { if (eta < -100) - rndEta = GetRandomEta(); + rndEta = GetRandomEta(kTRUE); if (phi < 0) - rndPhi = GetRandomPhi(); + rndPhi = GetRandomPhi(kTRUE); fGeom->GetAbsCellIdFromEtaPhi(rndEta, rndPhi, absId); repeats++; } while (absId == -1 && repeats < 100); @@ -727,19 +899,47 @@ void AliJetModelBaseTask::GetRandomCell(Double_t &eta, Double_t &phi, Int_t &abs } //________________________________________________________________________ -Double_t AliJetModelBaseTask::GetRandomEta() +Double_t AliJetModelBaseTask::GetRandomEta(Bool_t emcal) { // Get random eta. - return gRandom->Rndm() * (fEtaMax - fEtaMin) + fEtaMin; + Double_t etamax = fEtaMax; + Double_t etamin = fEtaMin; + + if (emcal) { + const Double_t EmcalMinEta = fGeom->GetArm1EtaMin(); + const Double_t EmcalMaxEta = fGeom->GetArm1EtaMax(); + + if (etamax > EmcalMaxEta) etamax = EmcalMaxEta; + if (etamax < EmcalMinEta) etamax = EmcalMinEta; + if (etamin > EmcalMaxEta) etamin = EmcalMaxEta; + if (etamin < EmcalMinEta) etamin = EmcalMinEta; + } + + return gRandom->Rndm() * (etamax - etamin) + etamin; } //________________________________________________________________________ -Double_t AliJetModelBaseTask::GetRandomPhi() +Double_t AliJetModelBaseTask::GetRandomPhi(Bool_t emcal) { // Get random phi. + + Double_t phimax = fPhiMax; + Double_t phimin = fPhiMin; + + if (emcal) { + const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad(); + const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad(); + + if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi; + if (phimax < EmcalMinPhi) phimax = EmcalMinPhi; + if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi; + if (phimin < EmcalMinPhi) phimin = EmcalMinPhi; + } - return gRandom->Rndm() * (fPhiMax - fPhiMin) + fPhiMin; + Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin; + + return result; } //________________________________________________________________________ @@ -753,6 +953,46 @@ Double_t AliJetModelBaseTask::GetRandomPt() return gRandom->Rndm() * (fPtMax - fPtMin) + fPtMin; } +//________________________________________________________________________ +void AliJetModelBaseTask::GetRandomParticle(Double_t &pt, Double_t &eta, Double_t &phi, Bool_t emcal) +{ + // Get a random particle. + + eta = GetRandomEta(emcal); + + if (fPtPhiEvPlDistribution) { + Double_t phimax = fPhiMax; + Double_t phimin = fPhiMin; + + if (emcal) { + const Double_t EmcalMinPhi = fGeom->GetArm1PhiMin() * TMath::DegToRad(); + const Double_t EmcalMaxPhi = fGeom->GetArm1PhiMax() * TMath::DegToRad(); + + if (phimax > EmcalMaxPhi) phimax = EmcalMaxPhi; + if (phimax < EmcalMinPhi) phimax = EmcalMinPhi; + if (phimin > EmcalMaxPhi) phimin = EmcalMaxPhi; + if (phimin < EmcalMinPhi) phimin = EmcalMinPhi; + } + + if (fPtPhiEvPlDistribution->GetXmin() > phimax || fPtPhiEvPlDistribution->GetXmax() < phimin) { + AliWarning(Form("The hisogram %s does not overlap with the EMCal acceptance limits. It will be ignored.",fPtPhiEvPlDistribution->GetName())); + pt = GetRandomPt(); + phi = GetRandomPhi(emcal); + } + else { + do { + fPtPhiEvPlDistribution->GetRandom2(pt,phi); + phi += fPsi; + if (phi > TMath::Pi() * 2) phi -= TMath::Pi() * 2; + } while (phi > phimax || phi < phimin); + } + } + else { + pt = GetRandomPt(); + phi = GetRandomPhi(emcal); + } +} + //________________________________________________________________________ void AliJetModelBaseTask::Run() {