1 // Simple class to generate toy events which can be fed into the jet finder
3 // Author: Redmer Alexander Bertens, Utrecht University, 2013
4 // rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
9 #include "TClonesArray.h"
12 #include "TParticle.h"
13 #include "TVirtualMCDecayer.h"
15 #include "AliAODEvent.h"
16 #include "AliAnalysisManager.h"
17 #include "AliAnalysisTaskJetFlowMC.h"
19 #include "AliPicoTrack.h"
21 class AliAnalysisTaskJetFlowMC;
24 ClassImp(AliAnalysisTaskJetFlowMC)
26 //_____________________________________________________________________________
27 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) {
28 // default constructor for root IO
29 for(Int_t i(0); i < 10; i++) {
34 fHistOriginalSpectrum[i] = 0x0;
35 fHistOriginalEtaPhi[i] = 0x0;
36 fHistToySpectrum[i] = 0x0;
37 fHistToyEtaPhi[i] = 0x0;
38 fHistOriginalDeltaPhi[i] = 0x0;
39 fHistToyDeltaPhi[i] = 0x0;
42 for(Int_t i(0); i < 25; i++) fDecayerCache[i] = 0x0;
44 //_____________________________________________________________________________
45 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) {
47 DefineInput(0, TChain::Class());
48 if(fQA) DefineOutput(1, TList::Class());
49 for(Int_t i(0); i < 10; i++) {
54 fHistOriginalSpectrum[i] = 0x0;
55 fHistOriginalEtaPhi[i] = 0x0;
56 fHistToySpectrum[i] = 0x0;
57 fHistToyEtaPhi[i] = 0x0;
58 fHistOriginalDeltaPhi[i] = 0x0;
59 fHistToyDeltaPhi[i] = 0x0;
62 // by default construction replace gRandom
64 if(gRandom) delete gRandom;
65 gRandom = new TRandom3(seed);
67 for(Int_t i(0); i < 25; i++) fDecayerCache[i] = 0x0;
69 //_____________________________________________________________________________
70 AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()
72 // desctructor, claim ownership of deleted objects by setting member pointers to NULL
80 for(Int_t i(0); i < 25; i++) {
81 if(fDecayerCache[i]) {
82 fDecayerCache[i]->Delete();
83 delete fDecayerCache[i];
84 fDecayerCache[i] = NULL;
88 fDecayerResults->Delete();
89 delete fDecayerResults;
90 fDecayerResults = NULL;
94 //_____________________________________________________________________________
95 void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()
97 // Create my user objects.
98 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
99 fTracksOut = new TClonesArray("AliPicoTrack");
100 fTracksOut->SetName(fTracksOutName);
102 fOutputList = new TList();
103 fOutputList->SetOwner(kTRUE);
105 if(!fCentralityClasses) { // classes must be defined at this point
106 Int_t c[] = {-10, 110};
107 fCentralityClasses = new TArrayI(sizeof(c)/sizeof(c[0]), c);
109 if(fHistIntV2 && !fHistIntV3) { // define function
110 fFuncVn = new TF1("kV2", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
111 fFuncVn->SetParameter(0, 1.); // normalization
112 fFuncVn->SetParameter(3, 0.2); // v2
113 fFuncVn->FixParameter(1, 1.); // constant
114 fFuncVn->FixParameter(2, 2.); // constant
115 } else if (!fHistIntV2 && fHistIntV3) {
116 fFuncVn = new TF1("kV3", "[0]*([1]+[2]*[3]*TMath::Cos([2]*(x-[4])))", 0, TMath::TwoPi());
117 fFuncVn->SetParameter(0, 1.); // normalization
118 fFuncVn->SetParameter(3, 0.2); // v3
119 fFuncVn->FixParameter(1, 1.); // constant
120 fFuncVn->FixParameter(2, 3.); // constant
121 } else if (fHistIntV2 && fHistIntV3) {
122 fFuncVn = new TF1("kCombined", "[0]*([1]+[2]*([3]*TMath::Cos([2]*(x-[4]))+[7]*TMath::Cos([5]*(x-[6]))))", 0, TMath::TwoPi());
123 fFuncVn->SetParameter(0, 1.); // normalization
124 fFuncVn->SetParameter(3, 0.2); // v2
125 fFuncVn->FixParameter(1, 1.); // constant
126 fFuncVn->FixParameter(2, 2.); // constant
127 fFuncVn->FixParameter(5, 3.); // constant
128 fFuncVn->SetParameter(7, 0.2); // v3
130 // add the generator objects that have been added to the task
131 if(fTrackSpectrum && fQA) fOutputList->Add(fTrackSpectrum);
132 if(fJetSpectrumSF && fQA) fOutputList->Add(fJetSpectrumSF);
133 if(fHistIntV2 && fQA) fOutputList->Add(fHistIntV2);
134 if(fHistIntV3 && fQA) fOutputList->Add(fHistIntV3);
135 // create the QA histos
137 for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
138 fHistOriginalSpectrum[i] = BookTH1F("fHistOriginalSpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
139 fHistOriginalEtaPhi[i] = BookTH2F("fHistOriginalEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
140 fHistToySpectrum[i] = BookTH1F("fHistToySpectrum", "p_{t} [GeV/c]", 200, 0, 200, i);
141 fHistToyEtaPhi[i] = BookTH2F("fHistToyEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi(), i);
142 fHistOriginalDeltaPhi[i] = BookTH1F("fHistOriginalDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
143 fHistToyDeltaPhi[i] = BookTH1F("fHistToyDeltaPhi", "#varphi - #Psi", 100, 0., TMath::Pi(), i);
144 fHistToyVn[i] = BookTH2F("fHistToyVn", "p_{t} [GeV/c]", "v_{n}", 1000, 0, 200, 200, 0, .2, i);
146 if(fFuncDiffV2[i]) fOutputList->Add(fFuncDiffV2[i]);
147 if(fFuncDiffV3[i]) fOutputList->Add(fFuncDiffV3[i]);
148 if(fHistDiffV2[i]) fOutputList->Add(fHistDiffV2[i]);
149 if(fHistDiffV3[i]) fOutputList->Add(fHistDiffV3[i]);
152 if(fJetSpectrumSF && fQA) {
153 fHistSFJetSpectrum = BookTH1F("fHistSFJetSpectrum", "p_{t} SF jets [GeV/c]", 100, 0, 200);
154 fHistSFJetEtaPhi = BookTH2F("fHistSFJetEtaPhi", "#eta", "#varphi", 100, -1., 1., 100, 0, TMath::TwoPi());
156 // reset the random seed
159 gRandom = new TRandom3(0);
162 // if presente initialize the decayer
165 // this may seem a bit amtibious but will reduce new / delete calls significantly
166 // to avoid memory fragmentation
167 for(Int_t i(0); i < 25; i++) fDecayerCache[i] = new TClonesArray("TParticle", 10);
168 if(fDecayerIterations > 3) {
169 printf(" -> Warning, fDecayerIterations set to %i. This number is too high, reverting to 3\n", fDecayerIterations);
170 fDecayerIterations = 3;
172 fDecayerResults = new TClonesArray("AliPicoTrack", 100);
176 PostData(1, fOutputList);
178 //_____________________________________________________________________________
179 TH1F* AliAnalysisTaskJetFlowMC::BookTH1F(const char* name, const char* x, Int_t bins, Double_t min, Double_t max, Int_t c, Bool_t append)
181 // book a TH1F and connect it to the output container
182 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
183 if(!fOutputList) return 0x0;
185 if(c!=-1) { // format centrality dependent histograms accordingly
186 name = Form("%s_%i", name, c);
187 title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
189 title += Form(";%s;[counts]", x);
190 TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
192 if(append) fOutputList->Add(histogram);
195 //_____________________________________________________________________________
196 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)
198 // book a TH2F and connect it to the output container
199 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
200 if(!fOutputList) return 0x0;
202 if(c!=-1) { // format centrality dependent histograms accordingly
203 name = Form("%s_%i", name, c);
204 title += Form("_%i-%i", fCentralityClasses->At(c), fCentralityClasses->At(1+c));
206 title += Form(";%s;%s", x, y);
207 TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
209 if(append) fOutputList->Add(histogram);
212 //_____________________________________________________________________________
213 void AliAnalysisTaskJetFlowMC::UserExec(Option_t *)
215 // user exec, called for each event.
216 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
217 if(!AliAnalysisManager::GetAnalysisManager()) return;
218 // retrieve tracks from input, only necessary when 'reuse' option is set, otherwise tracks will ge generated
219 if (fReuseTracks && !fTracksIn) {
220 fTracksIn = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject(fTracksInName));
222 AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data()));
226 // get the centrality bin for qa plots
227 Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
229 for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
230 if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {
234 if(fCenBin < 0) return;
235 // add tracks to event if not yet there
236 fTracksOut->Delete();
237 if (!(InputEvent()->FindListObject(fTracksOutName))) InputEvent()->AddObject(fTracksOut);
239 // get the event plane
240 CalculateEventPlane();
243 const Int_t Ntracks = fTracksIn->GetEntriesFast();
244 for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
245 AliPicoTrack* track = static_cast<AliPicoTrack*>(fTracksIn->At(iTracks));
247 Double_t phi(track->Phi()), pt((fTrackSpectrum) ? GetTrackPt() : track->Pt()), eta(fRandomizeEta ? GetTrackEta() : track->Eta());
248 // fill qa histo's before applying any (possible) afterburner
249 if(fQA) FillHistogramsOriginalData(pt, eta, phi);
250 if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);
251 else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);
252 else if(fHistIntV2 || fHistIntV3) SampleVnFromTF1(phi);
253 if(fDecayer) nacc = InsertDecayDaughters(track, fTracksOut);
255 /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1);
260 Double_t pt(0), eta(0), phi(0);
262 for (Int_t iTracks = 0; iTracks < fMult; ++iTracks) {
264 eta = gRandom->Uniform(-.9, .9);
265 phi = gRandom->Uniform(0., TMath::TwoPi());
266 charge = (gRandom->Uniform() < 0.5) ? -1 : 1;
267 // fill qa histo's before applying any (possible) afterburner
268 if(fQA) FillHistogramsOriginalData(pt, eta, phi);
269 if(fHistDiffV2[fCenBin] || fFuncDiffV2[fCenBin]) V2AfterBurner(phi, eta, pt);
270 else if(fHistDiffV3[fCenBin] || fFuncDiffV3[fCenBin]) V3AfterBurner(phi, eta, pt);
271 else if(fHistIntV2 || fHistIntV3) SampleVnFromTF1(phi);
272 if(fDecayer) nacc = InsertDecayDaughters(pt, phi, eta, -999, charge, fTracksOut);
274 /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, 1, 0, 4, eta, phi, pt, 1);
279 if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);
280 if(fQA) PostData(1, fOutputList);
281 if(fDebug > 0) PrintInfo();
283 //_____________________________________________________________________________
284 void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
286 // similar to AliFlowTrackSimple::AddV2, except for the flow fluctuations
287 if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
288 phi = gRandom->Uniform(0, TMath::TwoPi());
289 Double_t phi0(phi), v2(GetV2(pt)), f(0.), fp(0.), phiprev(0.);
290 if(TMath::AreEqualAbs(v2, 0, 1e-5) && fQA) {
291 FillHistogramsToyData(pt, eta, phi, v2);
294 // introduce flow fluctuations (gaussian)
295 if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
296 for (Int_t i=0; i!=fMaxNumberOfIterations; ++i) {
297 phiprev=phi; //store last value for comparison
298 f = phi-phi0+v2*TMath::Sin(2.*(phi-fPsi2));
299 fp = 1.0+2.0*v2*TMath::Cos(2.*(phi-fPsi2)); //first derivative
301 if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
303 if(fQA) FillHistogramsToyData(pt, eta, phi, v2);
305 //_____________________________________________________________________________
306 void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
308 // similar to AliFlowTrackSimple::AddV3, except for the flow fluctuations
309 if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
310 phi = gRandom->Uniform(0, TMath::TwoPi());
311 Double_t phi0(phi), v3(GetV3(pt)), f(0.), fp(0.), phiprev(0.);
312 if(TMath::AreEqualAbs(v3, 0, 1e-5) && fQA) {
313 FillHistogramsToyData(pt, eta, phi, v3);
316 // introduce flow fluctuations (gaussian)
317 if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
318 for (Int_t i=0; i<fMaxNumberOfIterations; i++) {
319 phiprev=phi; //store last value for comparison
320 f = phi-phi0+2./3.*v3*TMath::Sin(3.*(phi-fPsi3));
321 fp = 1.0+2.0*v3*TMath::Cos(3.*(phi-fPsi3)); //first derivative
323 if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
325 if(fQA) FillHistogramsToyData(pt, eta, phi, v3);
327 //_____________________________________________________________________________
328 void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)
330 // inject single fragmentation jet spectrum to the tclones array, note that emcal params
331 // equal the barrel kinematics to pass the track and jet cuts later on
332 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
333 for(Int_t i(nacc); i < (nacc + fNoOfSFJets); i++) {
334 Double_t eta(gRandom->Uniform(-.5, .5)), phi(gRandom->Uniform(0, TMath::TwoPi())), pt(fJetSpectrumSF->GetRandom());
335 /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[i]) AliPicoTrack(pt, eta, phi, +1, 0, 0, eta, phi, pt, 0);
337 fHistSFJetSpectrum->Fill(pt);
338 fHistSFJetEtaPhi->Fill(eta, phi);
344 //_____________________________________________________________________________
345 void AliAnalysisTaskJetFlowMC::CalculateEventPlane() {
346 // grab the event plane orientation from the AliVEvent header
347 if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
348 Double_t a(0), b(0), e(0), f(0);
349 switch (fDetectorType) {
351 fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);
352 fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);
356 fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);
357 fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);
361 fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;
362 fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);
366 // for fixed EP the EP is fixed to a constant value
373 // if requested, pass the event plane to the phi distribution
374 if(fHistIntV2 && fHistIntV3) {
375 fFuncVn->SetParameter(4, fPsi2); fFuncVn->SetParameter(6, fPsi3);
376 Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
377 Double_t v3(fHistIntV2->GetBinContent(fHistIntV3->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
378 if(fFlowFluctuations > -10) {
379 GetFlowFluctuation(v2); GetFlowFluctuation(v3);
381 fFuncVn->SetParameter(3, v2); fFuncVn->SetParameter(7, v3);
382 } else if (fHistIntV2) {
383 fFuncVn->SetParameter(4, fPsi2);
384 Double_t v2(fHistIntV2->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
385 if(fFlowFluctuations > -10) GetFlowFluctuation(v2);
386 fFuncVn->SetParameter(3, v2);
387 } else if (fHistIntV3) {
388 fFuncVn->SetParameter(4, fPsi3);
389 Double_t v3(fHistIntV3->GetBinContent(fHistIntV2->GetXaxis()->FindBin(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"))));
390 if(fFlowFluctuations > -10) GetFlowFluctuation(v3);
391 fFuncVn->SetParameter(3, v3);
394 //_____________________________________________________________________________
395 Int_t AliAnalysisTaskJetFlowMC::InsertDecayDaughters(
396 AliPicoTrack* mother,
397 TClonesArray* daughters)
399 // wrapper function to return decay daughters
400 return InsertDecayDaughters(
408 //_____________________________________________________________________________
409 Int_t AliAnalysisTaskJetFlowMC::InsertDecayDaughters(
415 TClonesArray* daughters)
417 // use TVirtualMCDecayer to force decay, daughters will be
418 // temporarily be stored in the TClonesArray as TParticle's
419 // but finally appended to the daughters TClonesArray array as AliPicoTracks
420 Int_t cacheCounter(0);
421 // bookkeeping variable, saves the status of the tclonesarray that is used
423 // [i] = correspond's to the i-th decay loop
424 // [fDecayerIterations] corresponds to the array that
425 // will eventually be returned (so all final state daughters)
426 static Int_t arrayStatus[25];
427 for(Int_t i(0); i < 25; i++) arrayStatus[i] = 0;
428 // allow for 'cascaded' decay, i.e. do multiple iterations
429 for(Int_t nestLevel(0); nestLevel < fDecayerIterations; nestLevel++) {
430 // for the first loop there is no TClonesArray available so grab the
431 // function arguments
433 TLorentzVector pMother( // 4-vector for mother
434 pt*TMath::Cos(phi), // px
435 pt*TMath::Sin(phi), // py
436 pt*TMath::SinH(eta), // pz
437 TMath::Sqrt(mass*mass+pt*pt*TMath::CosH(eta)*TMath::CosH(eta))); // total energy
438 // make sure input array is empty
439 fDecayerCache[cacheCounter]->Delete();
440 // FIXME we need a fix for the PDG code ...
441 Int_t pdgCode(211); // default to charged pion
442 if(mass == 0.13957) pdgCode = (charge > 0 ) ? 211 : -211;
443 // decay the mother and import the daughters
444 fDecayer->Decay(pdgCode, &pMother);
445 // fill first cache level with imported particles and set the decayer flag
446 fDecayer->ImportParticles(fDecayerCache[cacheCounter]);
447 arrayStatus[cacheCounter] = 1; // means that arrayStatus[0] is ready for decay in the next loop
451 // in subsequent loops check which arrays need to be decayed
452 for(Int_t j(0); j < 25; j++) {
453 // check if the array is supposed to be decayed at this nest level
454 if(arrayStatus[j] == nestLevel) {
455 // grab each particle from the array
456 for(Int_t k(0); k < fDecayerCache[j]->GetEntries(); k++) {
457 // 0 is the mother. only decay the mother if there are no daughters
458 if(k == 0 && fDecayerCache[j]->GetEntries() > 1) continue;
459 TParticle* daughter = static_cast<TParticle*>(fDecayerCache[j]->At(k));
460 if(!daughter) continue;
461 TLorentzVector pMother( // 4-vector for mother
465 daughter->Energy()); // total energy
466 // make sure input array is empty
467 fDecayerCache[cacheCounter]->Delete();
468 // FIXME we need a fix for the PDG code ...
469 // decay the mother and import the daughters
470 fDecayer->Decay(daughter->GetPdgCode(), &pMother);
471 fDecayer->ImportParticles(fDecayerCache[cacheCounter]);
472 // update the cache value
476 // update the next level
481 // what is left now is appending the freshly created tracks to the final track array
482 Int_t offset(daughters->GetEntries());
483 for(Int_t i(0); i < 25; i++) {
484 // only take the final state arrays
485 if(arrayStatus[i] == fDecayerIterations-1) {
486 for(Int_t j(0); j < fDecayerCache[j]->GetEntries(); j++) {
487 // if a mother has daughters, only push the daughters
488 if(j == 0 && fDecayerCache[j]->GetEntries() > 1) continue;
489 TParticle* daughter = static_cast<TParticle*>(fDecayerCache[j]->At(j));
490 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);
497 //_____________________________________________________________________________
498 void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)
502 //_____________________________________________________________________________
503 void AliAnalysisTaskJetFlowMC::PrintInfo() const
506 printf(" > No of retrieved tracks from %s \n \t %i \n", fTracksInName.Data(), fTracksIn->GetEntriesFast());
507 printf(" > No of created tracks in %s \n \t %i \n", fTracksOutName.Data(), fTracksOut->GetEntriesFast());
510 //_____________________________________________________________________________