#include <TLorentzVector.h>
#include <TRandom3.h>
#include <TList.h>
+#include <TF2.h>
#include "AliVEvent.h"
#include "AliAODCaloCluster.h"
fOutCellsName(),
fMCParticlesName(),
fOutMCParticlesName(),
+ fPythiaInfoName(""),
fIsMC(kFALSE),
fSuffix(),
fEtaMin(-1),
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),
fOutMCParticlesMap(0),
fMCLabelShift(0),
fEsdMode(kFALSE),
- fOutput(0)
+ fOutput(0),
+ fPythiaInfo(0x0)
{
// Default constructor.
+
+ fVertex[0] = 0;
+ fVertex[1] = 0;
+ fVertex[2] = 0;
}
//________________________________________________________________________
fOutCellsName(""),
fMCParticlesName(""),
fOutMCParticlesName(""),
+ fPythiaInfoName(""),
fIsMC(kFALSE),
fSuffix("Processed"),
fEtaMin(-1),
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),
fOutMCParticlesMap(0),
fMCLabelShift(0),
fEsdMode(kFALSE),
- fOutput(0)
+ fOutput(0),
+ fPythiaInfo(0x0)
{
// Standard constructor.
if (fQAhistos) {
DefineOutput(1, TList::Class());
}
+
+ fVertex[0] = 0;
+ fVertex[1] = 0;
+ fVertex[2] = 0;
}
//________________________________________________________________________
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();
if (fOutMCParticles)
fOutMCParticles->Delete();
}
+
+ if (fDensitySpectrum) {
+ fNTracks = TMath::Nint(fDensitySpectrum->GetRandom());
+ fNCells = TMath::Nint(fDensitySpectrum->GetRandom());
+ fNClusters = TMath::Nint(fDensitySpectrum->GetRandom());
+ }
// Clear map
if (fOutMCParticlesMap)
}
}
+ if (fPtPhiEvPlDistribution || fAddV2)
+ fPsi = gRandom->Rndm() * TMath::Pi();
+
Run();
-
- if (fCaloCells) {
+
+ if (fCaloCells && !fCopyArray) {
delete fCaloCells;
fCaloCells = tempCaloCells;
}
}
}
- if (fNTracks > 0 && !fTracksName.IsNull()) {
+ if (!fTracksName.IsNull()) {
fTracks = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksName));
if (!fTracks) {
AliWarning(Form("%s: Couldn't retrieve tracks with name %s!", GetName(), fTracksName.Data()));
}
}
- if (fNClusters > 0 && !fCaloName.IsNull()) {
+ if(fAddV2 && (!fDifferentialV2)) {
+ AliWarning(Form("%s: Cannot add v2 without diffential v2!", GetName()));
+ }
+
+ if (!fCaloName.IsNull()) {
fClusters = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fCaloName));
if (!fClusters) {
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 AliNamedArrayI(fMCParticlesName + "_Map", 9999);
- for (Int_t i = 0; i < 9999; i++) {
+ fMCParticlesMap = new AliNamedArrayI(fMCParticlesName + "_Map", 99999);
+ for (Int_t i = 0; i < 99999; i++) {
fMCParticlesMap->AddAt(i,i);
}
}
InputEvent()->AddObject(fOutMCParticles);
}
- fOutMCParticlesMap = new AliNamedArrayI(fOutMCParticlesName + "_Map",9999);
+ 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;
{
// Add a cell to the event.
- if (label == 0)
- label = fMarkMC;
- else
- label += fMCLabelShift;
+ if (label < 0)
+ label = 0;
+
+ label += fMarkMC + fMCLabelShift;
+
+ Short_t pos = -1;
+ if (fCaloCells)
+ pos = fCaloCells->GetCellPosition(absId);
- Short_t pos = fOutCaloCells->GetCellPosition(absId);
Double_t efrac = 1;
Bool_t increaseOnSuccess = kFALSE;
Double_t old_efrac = 0;
fOutCaloCells->GetCell(pos, cellNumber, old_e, old_time, old_label, old_efrac);
- if (e / (old_e + e) < 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;
}
- efrac = (e + old_e * old_efrac) / (e + old_e);
e += old_e;
}
//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 (nlabels != 0 && labels) {
- AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
- if (esdClus) {
- TArrayI parents(nlabels, labels);
- esdClus->AddLabels(parents);
- }
- else {
- AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
- if (aodClus)
- aodClus->SetLabel(labels, nlabels);
+ if (fMarkMC+fMCLabelShift != 0) {
+ for (UInt_t i = 0; i < nlabels; i++) {
+ parents[i] += fMarkMC+fMCLabelShift;
}
}
+ AliESDCaloCluster *esdClus = dynamic_cast<AliESDCaloCluster*>(dc);
+ if (esdClus) {
+ esdClus->AddLabels(parents);
+ }
+ else {
+ AliAODCaloCluster *aodClus = dynamic_cast<AliAODCaloCluster*>(dc);
+ if (aodClus)
+ aodClus->SetLabel(parents.GetArray(), nlabels);
+ }
+
return dc;
}
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<AliESDCaloCluster*>(cluster);
}
//________________________________________________________________________
-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;
}
AliAODMCParticle *aodpart = new ((*fOutMCParticles)[nPart]) AliAODMCParticle(*part);
+ 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.;
mclabel = 0;
fOutCaloCells->SetCell(i, cellNum, amp, time, mclabel, efrac);
+ fAddedCells++;
}
- fAddedCells = fCaloCells->GetNumberOfCells();
-
AliDebug(2, Form("%d cells from the current event", fAddedCells));
}
nCopiedPart++;
}
- if (!fMCParticlesMap)
+ if (!fMCParticlesMap || !fOutMCParticlesMap)
return;
+ if (fOutMCParticlesMap->GetSize() < fMCParticlesMap->GetSize())
+ fOutMCParticlesMap->Set(fMCParticlesMap->GetSize() * 2);
+
for (Int_t i = 0; i < fMCParticlesMap->GetSize(); i++) {
fOutMCParticlesMap->AddAt(fMCParticlesMap->At(i), i);
if (fMCParticlesMap->At(i) >= 0)
if (phimin < EmcalMinPhi) phimin = EmcalMinPhi;
}
- return gRandom->Rndm() * (phimax - phimin) + phimin;
+ Double_t result = gRandom->Rndm() * (phimax - phimin) + phimin;
+
+ return result;
}
//________________________________________________________________________
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()
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff