PWG/FLOW/Tasks/AliAnalysisTaskJetFlowMC.cxx

   1 // Simple class to generate toy events which can be fed into the jet finder
   2 //
   3 // Author: Redmer Alexander Bertens, Utrecht University, 2013
   4 // rbertens@cern.ch, rbertens@nikhef.nl, r.a.bertens@uu.nl
   5
   6 // root includes
   7 #include "TChain.h"
   8 #include "TList.h"
   9 #include "TClonesArray.h"
  10 #include "TArrayI.h"
  11 #include "TRandom3.h"
  12 #include "TParticle.h"
  13 #include "TVirtualMCDecayer.h"
  14 // aliroot includes
  15 #include "AliAODEvent.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliAnalysisTaskJetFlowMC.h"
  18 #include "AliLog.h"
  19 #include "AliPicoTrack.h"
  20
  21 class AliAnalysisTaskJetFlowMC;
  22 using namespace std;
  23
  24 ClassImp(AliAnalysisTaskJetFlowMC)
  25
  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++) {
  30         fFuncDiffV2[i]                  = 0x0;
  31         fFuncDiffV3[i]                  = 0x0;
  32         fHistDiffV2[i]                  = 0x0;
  33         fHistDiffV3[i]                  = 0x0;
  34         fHistOriginalSpectrum[i]        = 0x0;
  35         fHistOriginalEtaPhi[i]          = 0x0;
  36         fHistToySpectrum[i]             = 0x0;
  37         fHistToyEtaPhi[i]               = 0x0;
  38         fHistOriginalDeltaPhi[i]        = 0x0;
  39         fHistToyDeltaPhi[i]             = 0x0;
  40         fHistToyVn[i]                   = 0x0;
  41     }
  42     for(Int_t i(0); i < 25; i++) fDecayerCache[i] = 0x0;
  43 }
  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) {
  46     // constructor
  47     DefineInput(0, TChain::Class());
  48     if(fQA) DefineOutput(1, TList::Class());
  49     for(Int_t i(0); i < 10; i++) {
  50         fFuncDiffV2[i]                  = 0x0;
  51         fFuncDiffV3[i]                  = 0x0;
  52         fHistDiffV2[i]                  = 0x0;
  53         fHistDiffV3[i]                  = 0x0;
  54         fHistOriginalSpectrum[i]        = 0x0;
  55         fHistOriginalEtaPhi[i]          = 0x0;
  56         fHistToySpectrum[i]             = 0x0;
  57         fHistToyEtaPhi[i]               = 0x0;
  58         fHistOriginalDeltaPhi[i]        = 0x0;
  59         fHistToyDeltaPhi[i]             = 0x0;
  60         fHistToyVn[i]                   = 0x0;
  61     }
  62     // by default construction replace gRandom
  63     if(seed > -1) {
  64         if(gRandom) delete gRandom;
  65         gRandom = new TRandom3(seed);
  66     }
  67     for(Int_t i(0); i < 25; i++) fDecayerCache[i] = 0x0;
  68 }
  69 //_____________________________________________________________________________
  70 AliAnalysisTaskJetFlowMC::~AliAnalysisTaskJetFlowMC()
  71 {
  72     // desctructor, claim ownership of deleted objects by setting member pointers to NULL
  73     if(fOutputList) {
  74       delete fOutputList;
  75       fOutputList= NULL;
  76     }
  77     if(fDecayer) {
  78         delete fDecayer;
  79         fDecayer = 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;
  85             }
  86         }
  87         if(fDecayerResults) {
  88             fDecayerResults->Delete();
  89             delete fDecayerResults;
  90             fDecayerResults = NULL;
  91         }
  92     }
  93 }
  94 //_____________________________________________________________________________
  95 void AliAnalysisTaskJetFlowMC::UserCreateOutputObjects()
  96 {
  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);
 101     if(fQA) {
 102         fOutputList = new TList();
 103         fOutputList->SetOwner(kTRUE);
 104     }
 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);
 108     }
 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
 129     }
 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
 136     if(fQA) {
 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);
 145             // add to outputlist
 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]);
 150         }
 151     }
 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());
 155     }
 156     // reset the random seed
 157     if(gRandom) {
 158         delete gRandom;
 159         gRandom = new TRandom3(0);
 160     }
 161     if(!fQA) return;
 162     // if presente initialize the decayer
 163     if(fDecayer) {
 164         fDecayer->Init();
 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;
 171         }
 172         fDecayerResults = new TClonesArray("AliPicoTrack", 100);
 173     }
 174     // post output
 175     fOutputList->Sort();
 176     PostData(1, fOutputList);
 177 }
 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)
 180 {
 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;
 184     TString title(name);
 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));
 188     }
 189     title += Form(";%s;[counts]", x);
 190     TH1F* histogram = new TH1F(name, title.Data(), bins, min, max);
 191     histogram->Sumw2();
 192     if(append) fOutputList->Add(histogram);
 193     return histogram;
 194 }
 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)
 197 {
 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;
 201     TString title(name);
 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));
 205     }
 206     title += Form(";%s;%s", x, y);
 207     TH2F* histogram = new TH2F(name, title.Data(), binsx, minx, maxx, binsy, miny, maxy);
 208     histogram->Sumw2();
 209     if(append) fOutputList->Add(histogram);
 210     return histogram;
 211 }
 212 //_____________________________________________________________________________
 213 void AliAnalysisTaskJetFlowMC::UserExec(Option_t *)
 214 {
 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));
 221         if (!fTracksIn) {
 222           AliError(Form("Could not retrieve tracks %s!", fTracksInName.Data()));
 223           return;
 224         }
 225     }
 226     // get the centrality bin for qa plots
 227     Double_t cent(InputEvent()->GetCentrality()->GetCentralityPercentile("V0M"));
 228     fCenBin = -1;
 229     for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
 230         if(cent >= fCentralityClasses->At(i) && cent <= fCentralityClasses->At(1+i)) {
 231             fCenBin = i;
 232             break; }
 233     }
 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);
 238     fTracksOut->Clear();
 239     // get the event plane
 240     CalculateEventPlane();
 241     Int_t nacc(0);
 242     if(fReuseTracks) {
 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));
 246             if(!track) continue;
 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);
 254             else {
 255                 /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, track->Charge(), track->GetLabel(), 4, track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), track->GetTrackPtOnEMCal(), 1);
 256                 nacc++;
 257             }
 258         }
 259     } else {
 260         Double_t pt(0), eta(0), phi(0);
 261         Short_t charge(0);
 262         for (Int_t iTracks = 0; iTracks < fMult; ++iTracks) {
 263             pt = GetTrackPt();
 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);
 273             else {
 274                 /*AliPicoTrack *picotrack =*/ new ((*fTracksOut)[nacc]) AliPicoTrack(pt, eta, phi, 1, 0, 4, eta, phi, pt, 1);
 275                 nacc++;
 276             }
 277         }
 278     }
 279     if(fJetSpectrumSF && fNoOfSFJets > 0) InjectSingleFragmentationJetSpectrum(nacc);
 280     if(fQA) PostData(1, fOutputList);
 281     if(fDebug > 0) PrintInfo();
 282 }
 283 //_____________________________________________________________________________
 284 void AliAnalysisTaskJetFlowMC::V2AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
 285 {
 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);
 292         return;
 293     }
 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
 300         phi -= f/fp;
 301         if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
 302     }
 303     if(fQA) FillHistogramsToyData(pt, eta, phi, v2);
 304 }
 305 //_____________________________________________________________________________
 306 void AliAnalysisTaskJetFlowMC::V3AfterBurner(Double_t &phi, Double_t &eta, Double_t &pt) const
 307 {
 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);
 314         return;
 315     }
 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
 322         phi -= f/fp;
 323         if (TMath::AreEqualAbs(phiprev,phi,fPrecisionPhi)) break;
 324     }
 325     if(fQA) FillHistogramsToyData(pt, eta, phi, v3);
 326 }
 327 //_____________________________________________________________________________
 328 void AliAnalysisTaskJetFlowMC::InjectSingleFragmentationJetSpectrum(Int_t nacc)
 329 {
 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);
 336         if(fQA) {
 337             fHistSFJetSpectrum->Fill(pt);
 338             fHistSFJetEtaPhi->Fill(eta, phi);
 339         }
 340         ++i;
 341     }
 342     nacc = 0;
 343 }
 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) {
 350         case kVZEROA : {
 351             fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, e, f);
 352             fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, a, b);
 353             break;
 354         }
 355         case kVZEROC : {
 356             fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, e, f);
 357             fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, a, b);
 358             break;
 359         }
 360         case kVZEROComb : {
 361             fPsi2 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, e, f) ;
 362             fPsi3 = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, a, b);
 363             break;
 364         }
 365         case kFixedEP : {
 366             // for fixed EP the EP is fixed to a constant value
 367             fPsi2 = 0.;
 368             fPsi3 = 1.;
 369             break;
 370         }
 371         default : break;
 372     }
 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);
 380         }
 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);
 392     }
 393 }
 394 //_____________________________________________________________________________
 395 Int_t AliAnalysisTaskJetFlowMC::InsertDecayDaughters(
 396         AliPicoTrack* mother,
 397         TClonesArray* daughters)
 398 {
 399     // wrapper function to return decay daughters
 400     return InsertDecayDaughters(
 401             mother->Pt(),
 402             mother->Phi(),
 403             mother->Eta(),
 404             mother->M(),
 405             mother->Charge(),
 406             daughters);
 407 }
 408 //_____________________________________________________________________________
 409 Int_t AliAnalysisTaskJetFlowMC::InsertDecayDaughters(
 410         Double_t pt,
 411         Double_t phi,
 412         Double_t eta,
 413         Double_t mass,
 414         Short_t charge,
 415         TClonesArray* daughters)
 416 {
 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
 422     // [0] = unused
 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
 432         if(nestLevel == 0) {
 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
 448             nestLevel++;
 449             cacheCounter++;
 450         } else {
 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
 462                             daughter->Px(),
 463                             daughter->Py(),
 464                             daughter->Pz(),
 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
 473                         cacheCounter++;
 474                     }
 475                 }
 476                 // update the next level
 477                 nestLevel++;
 478             }
 479         }
 480     }
 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);
 491                 offset++;
 492             }
 493         }
 494     }
 495     return offset;
 496 }
 497 //_____________________________________________________________________________
 498 void AliAnalysisTaskJetFlowMC::Terminate(Option_t *)
 499 {
 500     // terminate
 501 }
 502 //_____________________________________________________________________________
 503 void AliAnalysisTaskJetFlowMC::PrintInfo() const
 504 {
 505     // print info
 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());
 508
 509 }
 510 //_____________________________________________________________________________