+++ /dev/null
-// Simple class to generate toy events which can be fed into the jet finder
-//
-// Author: Redmer Alexander Bertens, Utrecht University, 2013
-// rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
-
-// root includes
-#include "TChain.h"
-#include "TList.h"
-#include "TClonesArray.h"
-#include "TArrayI.h"
-#include "TRandom3.h"
-#include "TParticle.h"
-#include "TVirtualMCDecayer.h"
-// aliroot includes
-#include "AliAODEvent.h"
-#include "AliAnalysisManager.h"
-#include "AliAnalysisTaskJetFlowMC.h"
-#include "AliLog.h"
-#include "AliPicoTrack.h"
-
-class AliAnalysisTaskJetFlowMC;
-using namespace std;
-
-ClassImp(AliAnalysisTaskJetFlowMC)
-
-//_____________________________________________________________________________
-AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC() : AliAnalysisTaskSE("AliAnalysisTaskJetFlowMC"), fQA(kFALSE), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kFixedEP), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0), fDecayer(0x0), fDecayerIterations(1), fDecayerResults(0x0) {
- // default constructor for root IO
- for(Int_t i(0); i < 10; i++) {
- fFuncDiffV2[i] = 0x0;
- fFuncDiffV3[i] = 0x0;
- fHistDiffV2[i] = 0x0;
- fHistDiffV3[i] = 0x0;
- fHistOriginalSpectrum[i] = 0x0;
- fHistOriginalEtaPhi[i] = 0x0;
- fHistToySpectrum[i] = 0x0;
- fHistToyEtaPhi[i] = 0x0;
- fHistOriginalDeltaPhi[i] = 0x0;
- fHistToyDeltaPhi[i] = 0x0;
- fHistToyVn[i] = 0x0;
- }
- for(Int_t i(0); i < 25; i++) fDecayerCache[i] = 0x0;
-}
-//_____________________________________________________________________________
-AliAnalysisTaskJetFlowMC::AliAnalysisTaskJetFlowMC(const char *name, Bool_t qa, Int_t seed) : AliAnalysisTaskSE(name), fQA(qa), fTracksOutName("JetFlowMC"), fTracksInName("PicoTrack"), fTracksIn(0), fTracksOut(0), fReuseTracks(kFALSE), fMult(2200), fCenBin(-1), fCentralityClasses(0), fFuncVn(0), fOutputList(0), fTrackSpectrum(0), fRandomizeEta(kTRUE), fJetSpectrumSF(0), fNoOfSFJets(0), fHistIntV2(0), fHistIntV3(0), fFlowFluctuations(-10), fMaxNumberOfIterations(100), fPsi2(-10), fPsi3(-10), fPrecisionPhi(1e-10), fDetectorType(kFixedEP), fHistSFJetSpectrum(0), fHistSFJetEtaPhi(0), fDecayer(0x0), fDecayerIterations(1), fDecayerResults(0x0) {
- // constructor
- DefineInput(0, TChain::Class());
- if(fQA) DefineOutput(1, TList::Class());
- for(Int_t i(0); i < 10; i++) {
- fFuncDiffV2[i] = 0x0;
- fFuncDiffV3[i] = 0x0;
- fHistDiffV2[i] = 0x0;
- fHistDiffV3[i] = 0x0;
- fHistOriginalSpectrum[i] = 0x0;
- fHistOriginalEtaPhi[i] = 0x0;
- fHistToySpectrum[i] = 0x0;
- fHistToyEtaPhi[i] = 0x0;
- fHistOriginalDeltaPhi[i] = 0x0;
- fHistToyDeltaPhi[i] = 0x0;
- fHistToyVn[i] = 0x0;
- }
- // by default construction replace gRandom
- if(seed > -1) {
- if(gRandom) delete gRandom;
- gRandom = new TRandom3(seed);
- }
- for(Int_t i(0); i < 25; i++) fDecayerCache[i] = 0x0;
-}
-//_____________________________________________________________________________
-AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()
-{
- // desctructor, claim ownership of deleted objects by setting member pointers to NULL
- if(fOutputList) {
- delete fOutputList;
- fOutputList= NULL;
- }
- if(fDecayer) {
- delete fDecayer;
- fDecayer = NULL;
- for(Int_t i(0); i < 25; i++) {
- if(fDecayerCache[i]) {
- fDecayerCache[i]->Delete();
- delete fDecayerCache[i];
- fDecayerCache[i] = NULL;
- }
- }
- if(fDecayerResults) {
- fDecayerResults->Delete();
- delete fDecayerResults;
- fDecayerResults = NULL;
- }
- }
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()
-{
- // Create my user objects.
- if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- fTracksOut = new TClonesArray("AliPicoTrack");
- fTracksOut->SetName(fTracksOutName);
- if(fQA) {
- fOutputList = new TList();
- fOutputList->SetOwner(kTRUE);
- }
- if(!fCentralityClasses) { // classes must be defined at this point
- Int_t c[] = {-10, 110};
- fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);
- }
- if(fHistIntV2 && !fHistIntV3) { // define function
- fFuncVn = new TF1("kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
- fFuncVn->SetParameter(0, 1.); // normalization
- fFuncVn->SetParameter(3, 0.2); // v2
- fFuncVn->FixParameter(1, 1.); // constant
- fFuncVn->FixParameter(2, 2.); // constant
- } else if (!fHistIntV2 && fHistIntV3) {
- fFuncVn = new TF1("kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
- fFuncVn->SetParameter(0, 1.); // normalization
- fFuncVn->SetParameter(3, 0.2); // v3
- fFuncVn->FixParameter(1, 1.); // constant
- fFuncVn->FixParameter(2, 3.); // constant
- } else if (fHistIntV2 && fHistIntV3) {
- fFuncVn = new TF1("kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());
- fFuncVn->SetParameter(0, 1.); // normalization
- fFuncVn->SetParameter(3, 0.2); // v2
- fFuncVn->FixParameter(1, 1.); // constant
- fFuncVn->FixParameter(2, 2.); // constant
- fFuncVn->FixParameter(5, 3.); // constant
- fFuncVn->SetParameter(7, 0.2); // v3
- }
- // add the generator objects that have been added to the task
- if(fTrackSpectrum && fQA) fOutputList->Add(fTrackSpectrum);
- if(fJetSpectrumSF && fQA) fOutputList->Add(fJetSpectrumSF);
- if(fHistIntV2 && fQA) fOutputList->Add(fHistIntV2);
- if(fHistIntV3 && fQA) fOutputList->Add(fHistIntV3);
- // create the QA histos
- if(fQA) {
- for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
- fHistOriginalSpectrum[i] = BookTH1F("fHistOriginalSpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
- fHistOriginalEtaPhi[i] = BookTH2F("fHistOriginalEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
- fHistToySpectrum[i] = BookTH1F("fHistToySpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
- fHistToyEtaPhi[i] = BookTH2F("fHistToyEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
- fHistOriginalDeltaPhi[i] = BookTH1F("fHistOriginalDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
- fHistToyDeltaPhi[i] = BookTH1F("fHistToyDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
- fHistToyVn[i] = BookTH2F("fHistToyVn", "p_{t} [GeV/c]", "v_{n}", 1000, 0, 200, 200, 0, .2, i);
- // add to outputlist
- if(fFuncDiffV2[i]) fOutputList->Add(fFuncDiffV2[i]);
- if(fFuncDiffV3[i]) fOutputList->Add(fFuncDiffV3[i]);
- if(fHistDiffV2[i]) fOutputList->Add(fHistDiffV2[i]);
- if(fHistDiffV3[i]) fOutputList->Add(fHistDiffV3[i]);
- }
- }
- if(fJetSpectrumSF && fQA) {
- fHistSFJetSpectrum = BookTH1F("fHistSFJetSpectrum", "p_{t} SF jets [GeV/c]", 100, 0, 200);
- fHistSFJetEtaPhi = BookTH2F("fHistSFJetEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi());
- }
- // reset the random seed
- if(gRandom) {
- delete gRandom;
- gRandom = new TRandom3(0);
- }
- if(!fQA) return;
- // if presente initialize the decayer
- if(fDecayer) {
- fDecayer->Init();
- // this may seem a bit amtibious but will reduce new / delete calls significantly
- // to avoid memory fragmentation
- for(Int_t i(0); i < 25; i++) fDecayerCache[i] = new TClonesArray("TParticle", 10);
- if(fDecayerIterations > 3) {
- printf(" -> Warning, fDecayerIterations set to %i. This number is too high, reverting to 3\n", fDecayerIterations);
- fDecayerIterations = 3;
- }
- fDecayerResults = new TClonesArray("AliPicoTrack", 100);
- }
- // post output
- fOutputList->Sort();
- PostData(1, fOutputList);
-}
-//_____________________________________________________________________________
-TH1F* AliAnalysisTaskJetFlowMC::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
-{
- // book a TH1F and connect it to the output container
- if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- if(!fOutputList) return 0x0;
- TString title(name);
- if(c!=-1) { // format centrality dependent histograms accordingly
- name = Form("%s_%i", name, c);
- title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
- }
- title += Form(";%s;[counts]", x);
- TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
- histogram->Sumw2();
- if(append) fOutputList->Add(histogram);
- return histogram;
-}
-//_____________________________________________________________________________
-TH2F* AliAnalysisTaskJetFlowMC::BookTH2F(const char* name, const char* x, const char*y, Int_t binsx, Double_t minx, Double_t maxx, Int_t binsy, Double_t miny, Double_t maxy, Int_t c, Bool_t append)
-{
- // book a TH2F and connect it to the output container
- if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- if(!fOutputList) return 0x0;
- TString title(name);
- if(c!=-1) { // format centrality dependent histograms accordingly
- name = Form("%s_%i", name, c);
- title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
- }
- title += Form(";%s;%s", x, y);
- TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
- histogram->Sumw2();
- if(append) fOutputList->Add(histogram);
- return histogram;
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::UserExec(Option_t *)
-{
- // user exec, called for each event.
- if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- if(!AliAnalysisManager::GetAnalysisManager()) return;
- // retrieve tracks from input, only necessary when 'reuse' option is set, otherwise tracks will ge generated
- if (fReuseTracks && !fTracksIn) {
- fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));
- if (!fTracksIn) {
- AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data()));
- return;
- }
- }
- // get the centrality bin for qa plots
- Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
- fCenBin = -1;
- for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
- if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {
- fCenBin = i;
- break; }
- }
- if(fCenBin < 0) return;
- // add tracks to event if not yet there
- fTracksOut->Delete();
- if (!(InputEvent()->FindListObject(fTracksOutName))) InputEvent()->AddObject(fTracksOut);
- fTracksOut->Clear();
- // get the event plane
- CalculateEventPlane();
- Int_t nacc(0);
- if(fReuseTracks) {
- const Int_t Ntracks = fTracksIn->GetEntriesFast();
- for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
- AliPicoTrack* track = static_cast<AliPicoTrack*>(fTracksIn->At(iTracks));
- if(!track) continue;
- Double_t phi(track->Phi()), pt((fTrackSpectrum) ? GetTrackPt() : track->Pt()), eta(fRandomizeEta ? GetTrackEta() : track->Eta());
- // fill qa histo's before applying any (possible) afterburner
- if(fQA) FillHistogramsOriginalData(pt, eta, phi);
- if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);
- else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);
- else if(fHistIntV2 || fHistIntV3) SampleVnFromTF1(phi);
- if(fDecayer) nacc = InsertDecayDaughters(track, fTracksOut);
- else {
- /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1);
- nacc++;
- }
- }
- } else {
- Double_t pt(0), eta(0), phi(0);
- Short_t charge(0);
- for (Int_t iTracks = 0; iTracks < fMult; ++iTracks) {
- pt = GetTrackPt();
- eta = gRandom->Uniform(-.9, .9);
- phi = gRandom->Uniform(0., TMath::TwoPi());
- charge = (gRandom->Uniform() < 0.5) ? -1 : 1;
- // fill qa histo's before applying any (possible) afterburner
- if(fQA) FillHistogramsOriginalData(pt, eta, phi);
- if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);
- else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);
- else if(fHistIntV2 || fHistIntV3) SampleVnFromTF1(phi);
- if(fDecayer) nacc = InsertDecayDaughters(pt, phi, eta, -999, charge, fTracksOut);
- else {
- /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, 1, 0, 4, eta, phi, pt, 1);
- nacc++;
- }
- }
- }
- if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);
- if(fQA) PostData(1, fOutputList);
- if(fDebug > 0) PrintInfo();
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
-{
- // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations
- if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- phi = gRandom->Uniform(0, TMath::TwoPi());
- Double_t phi0(phi), v2(GetV2(pt)), f(0.), fp(0.), phiprev(0.);
- if(TMath::AreEqualAbs(v2, 0, 1e-5) && fQA) {
- FillHistogramsToyData(pt, eta, phi, v2);
- return;
- }
- // introduce flow fluctuations (gaussian)
- if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
- for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {
- phiprev=phi; //store last value for comparison
- f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));
- fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative
- phi -= f/fp;
- if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
- }
- if(fQA) FillHistogramsToyData(pt, eta, phi, v2);
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
-{
- // similar to AliFlowTrackSimple::AddV3, except for the flow fluctuations
- if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- phi = gRandom->Uniform(0, TMath::TwoPi());
- Double_t phi0(phi), v3(GetV3(pt)), f(0.), fp(0.), phiprev(0.);
- if(TMath::AreEqualAbs(v3, 0, 1e-5) && fQA) {
- FillHistogramsToyData(pt, eta, phi, v3);
- return;
- }
- // introduce flow fluctuations (gaussian)
- if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
- for (Int_t i=0; i<fMaxNumberOfIterations; i++) {
- phiprev=phi; //store last value for comparison
- f = phi-phi0+2./3.*v3*TMath::Sin(3.*(phi-fPsi3));
- fp = 1.0+2.0*v3*TMath::Cos(3.*(phi-fPsi3)); //first derivative
- phi -= f/fp;
- if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
- }
- if(fQA) FillHistogramsToyData(pt, eta, phi, v3);
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)
-{
- // inject single fragmentation jet spectrum to the tclones array, note that emcal params
- // equal the barrel kinematics to pass the track and jet cuts later on
- if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- for(Int_t i(nacc); i < (nacc + fNoOfSFJets); i++) {
- Double_t eta(gRandom->Uniform(-.5, .5)), phi(gRandom->Uniform(0, TMath::TwoPi())), pt(fJetSpectrumSF->GetRandom());
- /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[i]) AliPicoTrack(pt, eta, phi, +1, 0, 0, eta, phi, pt, 0);
- if(fQA) {
- fHistSFJetSpectrum->Fill(pt);
- fHistSFJetEtaPhi->Fill(eta, phi);
- }
- ++i;
- }
- nacc = 0;
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::CalculateEventPlane() {
- // grab the event plane orientation from the AliVEvent header
- if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
- Double_t a(0), b(0), e(0), f(0);
- switch (fDetectorType) {
- case kVZEROA : {
- fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);
- fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);
- break;
- }
- case kVZEROC : {
- fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);
- fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);
- break;
- }
- case kVZEROComb : {
- fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;
- fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);
- break;
- }
- case kFixedEP : {
- // for fixed EP the EP is fixed to a constant value
- fPsi2 = 0.;
- fPsi3 = 1.;
- break;
- }
- default : break;
- }
- // if requested, pass the event plane to the phi distribution
- if(fHistIntV2 && fHistIntV3) {
- fFuncVn->SetParameter(4, fPsi2); fFuncVn->SetParameter(6, fPsi3);
- Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
- Double_t v3(fHistIntV2->GetBinContent(fHistIntV3->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
- if(fFlowFluctuations > -10) {
- GetFlowFluctuation(v2); GetFlowFluctuation(v3);
- }
- fFuncVn->SetParameter(3, v2); fFuncVn->SetParameter(7, v3);
- } else if (fHistIntV2) {
- fFuncVn->SetParameter(4, fPsi2);
- Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
- if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
- fFuncVn->SetParameter(3, v2);
- } else if (fHistIntV3) {
- fFuncVn->SetParameter(4, fPsi3);
- Double_t v3(fHistIntV3->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
- if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
- fFuncVn->SetParameter(3, v3);
- }
-}
-//_____________________________________________________________________________
-Int_t AliAnalysisTaskJetFlowMC::InsertDecayDaughters(
- AliPicoTrack* mother,
- TClonesArray* daughters)
-{
- // wrapper function to return decay daughters
- return InsertDecayDaughters(
- mother->Pt(),
- mother->Phi(),
- mother->Eta(),
- mother->M(),
- mother->Charge(),
- daughters);
-}
-//_____________________________________________________________________________
-Int_t AliAnalysisTaskJetFlowMC::InsertDecayDaughters(
- Double_t pt,
- Double_t phi,
- Double_t eta,
- Double_t mass,
- Short_t charge,
- TClonesArray* daughters)
-{
- // use TVirtualMCDecayer to force decay, daughters will be
- // temporarily be stored in the TClonesArray as TParticle's
- // but finally appended to the daughters TClonesArray array as AliPicoTracks
- Int_t cacheCounter(0);
- // bookkeeping variable, saves the status of the tclonesarray that is used
- // [0] = unused
- // [i] = correspond's to the i-th decay loop
- // [fDecayerIterations] corresponds to the array that
- // will eventually be returned (so all final state daughters)
- static Int_t arrayStatus[25];
- for(Int_t i(0); i < 25; i++) arrayStatus[i] = 0;
- // allow for 'cascaded' decay, i.e. do multiple iterations
- for(Int_t nestLevel(0); nestLevel < fDecayerIterations; nestLevel++) {
- // for the first loop there is no TClonesArray available so grab the
- // function arguments
- if(nestLevel == 0) {
- TLorentzVector pMother( // 4-vector for mother
- pt*TMath::Cos(phi), // px
- pt*TMath::Sin(phi), // py
- pt*TMath::SinH(eta), // pz
- TMath::Sqrt(mass*mass+pt*pt*TMath::CosH(eta)*TMath::CosH(eta))); // total energy
- // make sure input array is empty
- fDecayerCache[cacheCounter]->Delete();
- // FIXME we need a fix for the PDG code ...
- Int_t pdgCode(211); // default to charged pion
- if(mass == 0.13957) pdgCode = (charge > 0 ) ? 211 : -211;
- // decay the mother and import the daughters
- fDecayer->Decay(pdgCode, &pMother);
- // fill first cache level with imported particles and set the decayer flag
- fDecayer->ImportParticles(fDecayerCache[cacheCounter]);
- arrayStatus[cacheCounter] = 1; // means that arrayStatus[0] is ready for decay in the next loop
- nestLevel++;
- cacheCounter++;
- } else {
- // in subsequent loops check which arrays need to be decayed
- for(Int_t j(0); j < 25; j++) {
- // check if the array is supposed to be decayed at this nest level
- if(arrayStatus[j] == nestLevel) {
- // grab each particle from the array
- for(Int_t k(0); k < fDecayerCache[j]->GetEntries(); k++) {
- // 0 is the mother. only decay the mother if there are no daughters
- if(k == 0 && fDecayerCache[j]->GetEntries() > 1) continue;
- TParticle* daughter = static_cast<TParticle*>(fDecayerCache[j]->At(k));
- if(!daughter) continue;
- TLorentzVector pMother( // 4-vector for mother
- daughter->Px(),
- daughter->Py(),
- daughter->Pz(),
- daughter->Energy()); // total energy
- // make sure input array is empty
- fDecayerCache[cacheCounter]->Delete();
- // FIXME we need a fix for the PDG code ...
- // decay the mother and import the daughters
- fDecayer->Decay(daughter->GetPdgCode(), &pMother);
- fDecayer->ImportParticles(fDecayerCache[cacheCounter]);
- // update the cache value
- cacheCounter++;
- }
- }
- // update the next level
- nestLevel++;
- }
- }
- }
- // what is left now is appending the freshly created tracks to the final track array
- Int_t offset(daughters->GetEntries());
- for(Int_t i(0); i < 25; i++) {
- // only take the final state arrays
- if(arrayStatus[i] == fDecayerIterations-1) {
- for(Int_t j(0); j < fDecayerCache[j]->GetEntries(); j++) {
- // if a mother has daughters, only push the daughters
- if(j == 0 && fDecayerCache[j]->GetEntries() > 1) continue;
- TParticle* daughter = static_cast<TParticle*>(fDecayerCache[j]->At(j));
- if(daughter) /*AliPicoTrack *picotrack =*/ new ((*daughters)[offset]) AliPicoTrack(daughter->Pt(), daughter->Eta(), daughter->Phi(), (daughter->GetPdgCode() > 0) ? 1 : -1, 0, 0, daughter->Eta(), daughter->Phi(), daughter->Pt(), 0);
- offset++;
- }
- }
- }
- return offset;
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)
-{
- // terminate
-}
-//_____________________________________________________________________________
-void AliAnalysisTaskJetFlowMC::PrintInfo() const
-{
- // print info
- printf(" > No of retrieved tracks from %s \n \t %i \n", fTracksInName.Data(), fTracksIn->GetEntriesFast());
- printf(" > No of created tracks in %s \n \t %i \n", fTracksOutName.Data(), fTracksOut->GetEntriesFast());
-
-}
-//_____________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff