/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/*****************************************************************
  AliFlowOnTheFlyEventGenerator

  Create events with realstic spectra and v_n values
  store as AliFlowEventSimple
  
  usage: see $ALICE_ROOT/PWGCF/FLOW/macros/GenerateOnTheFlyEvents.C

  origin:       Redmer Alexander Bertens, rbertens@cern.ch, rbertens@nikhef.nl, rbertens@uu.nl
                Carlos Eugenio Perez Lara
                Andrea Dubla
******************************************************************/

#include "Riostream.h"
#include "TObjArray.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TList.h"
#include "TTree.h"
#include "TParticle.h"
#include "TMath.h"
#include "TH1.h"
#include "TH1F.h"
#include "TH1D.h"
#include "TF1.h"
#include "TString.h"
#include "TProfile.h"
#include "TParameter.h"
#include "TBrowser.h"
#include "AliFlowTrackSimple.h"
#include "AliFlowEventSimple.h"
#include "TRandom.h"
#include "TVirtualMCDecayer.h"
#include "AliFlowOnTheFlyEventGenerator.h"

using std::endl;
using std::cout;

static Int_t OnTheFlyEventGeneratorCounter = 0;

ClassImp(AliFlowOnTheFlyEventGenerator)
//_____________________________________________________________________________
AliFlowOnTheFlyEventGenerator::AliFlowOnTheFlyEventGenerator() : fGenerators(0), fEmbedMe(0), fFlowEvent(0), fDecayer(0), fAddV2Mothers(0), fAddV3Mothers(0), fAddV2Daughters(0), fAddV3Daughters(0), fPsi2(0), fPsi3(0), fPrecisionPhi(.001), fMaxNumberOfIterations(100), fQA(0), fFF(0) {/* constructor used for IO by root */}
//_____________________________________________________________________________
AliFlowOnTheFlyEventGenerator::AliFlowOnTheFlyEventGenerator(Bool_t qa, Int_t ff, Int_t mult, TVirtualMCDecayer* decayer, Bool_t a, Bool_t b, Bool_t c, Bool_t d) : fGenerators(0), fEmbedMe(0), fFlowEvent(0), fDecayer(0), fAddV2Mothers(0), fAddV3Mothers(0), fAddV2Daughters(0), fAddV3Daughters(0), fPsi2(0), fPsi3(0), fPrecisionPhi(.001), fMaxNumberOfIterations(100), fQA(qa), fFF(ff) {
    // contructor
    if(!fFlowEvent) fFlowEvent = new AliFlowEventSimple(mult, (AliFlowEventSimple::ConstructionMethod)0, 0x0, 0, TMath::TwoPi(), -.9, .9);
    if(!fGenerators) {
        fGenerators = new TObjArray();
        fGenerators->SetOwner(kTRUE);
        InitGenerators();
    }
    fDecayer = decayer;         // decayer: user has ownership of decayer (see dtor)
    if(fDecayer) fDecayer->Init();
    fAddV2Mothers = a;
    fAddV3Mothers = b;
    fAddV2Daughters = c; 
    fAddV3Daughters = d;
    if(fAddV3Mothers||fAddV3Daughters) printf(" > WARNING - v3 is not impelmented yet ! \n <");
}
//_____________________________________________________________________________
AliFlowOnTheFlyEventGenerator::~AliFlowOnTheFlyEventGenerator() 
{
    // destructor
    if(fGenerators)     delete fGenerators; 
    if(fFlowEvent)      delete fFlowEvent;
}
//_____________________________________________________________________________
AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator* AliFlowOnTheFlyEventGenerator::Find(Short_t pdg, Bool_t make)  
{ 
    // return instance of generator class for given fPDG value
    for(int i(0); i < fGenerators->GetEntriesFast(); i++) {
        if(((NaiveFlowAndSpectrumGenerator*)fGenerators->At(i))->GetPDGCode()==pdg) {
            return (NaiveFlowAndSpectrumGenerator*)(fGenerators->At(i));
        }
    }
    if(make) {
        AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator* _tmp = new AliFlowOnTheFlyEventGenerator::NaiveFlowAndSpectrumGenerator(pdg, fQA, fFF);
        fGenerators->Add(_tmp);
        return _tmp;
    }
    return 0x0;
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::SetPtSpectrum(const char* func, Short_t species)
{
    // set pt spectrum for species, override default
    NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
    TF1* spectrum = gen->GetPtSpectrum();
    TString name = "";
    if(spectrum) {
        name = spectrum->GetName();
        delete spectrum;
    }
    else name = Form("pt_%i", (int)species);
    gen->SetPtSpectrum(new TF1(name.Data(), func, 0., 20));
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::SetPtDependentV2(const char* func, Short_t species)
{
    // set pt dependent v2 for species, override default
    NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
    TF1* v2 = gen->GetDifferentialV2();
    TString name = "";
    if(v2) { // if v2 is already present, get its name and release its memory
        name = v2->GetName();
        delete v2;
    }
    else name = Form("v2_%i", (int)species);
    gen->SetDifferentialV2(new TF1(name.Data(),func, 0., 20));
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::SetPtDependentV3(const char* func, Short_t species)
{
    // set pt dependent v2 for species, override default
    NaiveFlowAndSpectrumGenerator* gen = Find(species, kTRUE);
    TF1* v3 = gen->GetDifferentialV3();
    TString name = "";
    if(v3) { // if v3 is already present, get its name and release its memory
        name = v3->GetName();
        delete v3;
    }
    else name = Form("v3_%i", (int)species);
    gen->SetDifferentialV3(new TF1(name.Data(), func, 0., 20));
}
//_____________________________________________________________________________
AliFlowEventSimple* AliFlowOnTheFlyEventGenerator::GenerateOnTheFlyEvent(TClonesArray* event, Int_t nSpecies, Int_t species[], Int_t mult[], Int_t bg, Bool_t fluc)
{
    // generate a flow event according to settings provided in GENERATOR INTERFACE
    fPsi2 = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
    // generate a random number which will be used for flow fuctuations (if requested)
    Double_t fluc_poi(-1), fluc_rp(-1);
    if(fFF == 1 || fFF == 3 ) fluc_poi = gRandom->Uniform();
    if(fFF == 2) fluc_rp = gRandom->Uniform();
    if(fFF == 3) fluc_rp = fluc_poi;    // poi and rp fluctuations are fully correlated
    // calculate the total multiplicity for the event
    Int_t multPiKPr[] = {(int)(.7*bg/2.), (int)(.2*bg/2.), (int)(.1*bg/2.)};
    Int_t fPDGPiKPr[] = {211, 321, 2212};
    Int_t totalRPMultiplicity(0), totalPOIMultiplicity(0);
    for(int i(0); i < nSpecies; i++) totalPOIMultiplicity+=mult[i];
    for(int i(0); i < 3; i++)        totalRPMultiplicity+=multPiKPr[i];
    Int_t totalMultiplicity(totalRPMultiplicity+totalPOIMultiplicity);
    // generate the particles of interest. if requested, a vn boost is given to the primary particles
    for(int i(0); i < nSpecies; i++) {
      if(fluc) GenerateOnTheFlyTracks((Int_t)(gRandom->Uniform(mult[i]-TMath::Sqrt(mult[i]), mult[i]+TMath::Sqrt(mult[i]))), species[i], event, fluc_poi);
        else GenerateOnTheFlyTracks(mult[i], species[i], event, fluc_poi);
    }
    // generate reference particles. if requested, a vn boost is given to the primary particles
    for(int i(0); i < 3; i++) {
        if(fluc) {
	  GenerateOnTheFlyTracks((Int_t)(gRandom->Uniform(multPiKPr[i]-TMath::Sqrt(multPiKPr[i]), multPiKPr[i]+TMath::Sqrt(mult[i]))), fPDGPiKPr[i], event, fluc_rp);
	  GenerateOnTheFlyTracks((Int_t)(gRandom->Uniform(multPiKPr[i]-TMath::Sqrt(multPiKPr[i]), multPiKPr[i]+TMath::Sqrt(mult[i]))), -1*fPDGPiKPr[i], event, fluc_rp);
        }
        else {
            GenerateOnTheFlyTracks(multPiKPr[i], fPDGPiKPr[i], event, fluc_rp);
            GenerateOnTheFlyTracks(multPiKPr[i], -1*fPDGPiKPr[i], event, fluc_rp);
        }
    }
    // if requested, decay the event
    if(fDecayer)        DecayOnTheFlyTracks(event); 
    // if an embedded event is given to the generator at this point, embed it
    if(fEmbedMe) {
        event->AbsorbObjects(fEmbedMe); // event has ownership of embedded event
        fEmbedMe = 0x0;                 // reset to aviod embeddding more than once
    }
    // if requested, the secondaries are given a vn boost
    if(fAddV2Daughters) AddV2(event);
    // convert the event to an aliflowsimple event
    // the tagging (rp) is done in this step, all tracks are tagged as rp's. 
    // daughters are made orphans (i.e. primaries with secondaries are rejected from the sample)
    // converting clears the event tclones array 
    return ConvertTClonesToFlowEvent(event, totalMultiplicity);
}
//_____________________________________________________________________________
AliFlowEventSimple* AliFlowOnTheFlyEventGenerator::ConvertTClonesToFlowEvent(TClonesArray* event, Int_t totalMultiplicity)
{
    // convert the content of a tclones array to an aliflowevent simple for output and tag particles as poi's and rp's
    // first step, see if there's already a flow event available in memory. if not create one
    if(!fFlowEvent) fFlowEvent = new AliFlowEventSimple(totalMultiplicity, (AliFlowEventSimple::ConstructionMethod)0, 0x0, 0, TMath::TwoPi(), -.9, .9);
    // if there is a flow event, clear the members without releasing the allocated memory
    else fFlowEvent->ClearFast();
    fFlowEvent->SetMCReactionPlaneAngle(fPsi2);
    // prepare a track
    TParticle* pParticle = 0x0;
    Int_t iSelParticlesRP(0);
    for (Int_t i(0); i < event->GetEntries(); i++) {    // begin the loop over the input tracks
        pParticle = (TParticle*)event->At(i);           // get the particle 
        if (!pParticle) continue;                       // skip if empty slot (no particle)
        if (pParticle->GetNDaughters()!=0) continue;    // see if the particle has daughters (if so, reject it)      
        AliFlowTrackSimple* pTrack = new AliFlowTrackSimple(pParticle);         // allocate space for the new flow track
        pTrack->SetWeight(pParticle->Pz());                                     // ugly hack: store pz here ...
        pTrack->SetID(pParticle->GetPdgCode());                                 // set pid code as id
        pTrack->SetForRPSelection(kTRUE);                                       // tag ALL particles as RP's, 
        // ugly hack 2: set charge to -1 for primaries and 1 for secondaries
        if(pParticle->GetFirstMother()==-1) pTrack->SetCharge(-1);
        else pTrack->SetCharge(1);
        iSelParticlesRP++;
        fFlowEvent->AddTrack(pTrack);
    }
    fFlowEvent->SetNumberOfRPs(iSelParticlesRP);
    // all trakcs have safely been copied so we can clear the event
    event->Clear();
    OnTheFlyEventGeneratorCounter = 0;
    return fFlowEvent;
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::AddV2(TParticle* part, Double_t v2, Double_t fluc)
{
    // afterburner, adds v2, uses Newton-Raphson iteration
    Double_t phi = part->Phi();
    Double_t phi0=phi;
    Double_t f=0.;
    Double_t fp=0.;
    Double_t phiprev=0.;
    // introduce flow fluctuations (gaussian)
    if(fluc > -.5) {
        v2+=TMath::Sqrt(2*(v2*.25)*(v2*.25))*TMath::ErfInverse(2*fluc-1);
        if(fQA) Find(part->GetPdgCode(), kTRUE)->FillV2(part->Pt(), 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;
    }
    part->SetMomentum( part->Pt()*TMath::Cos(phi), part->Pt()*TMath::Sin(phi), part->Pz(), part->Energy() );
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::AddV2(TClonesArray* event)
{
    // afterburner, adds v2 for different species to all tracks in an event
    Double_t fluc(gRandom->Uniform());  // get a random number in case of fluctuations
    // FIXME at the moment no distincition between mothers and daughters
    TParticle *part;
    for(Int_t nTrack=0; nTrack!=event->GetEntriesFast(); ++nTrack) {
        part = (TParticle*) event->At(nTrack);
        // for each track, call overloaded addv2 function with the correct generator
        // create a generator in the case where the decayer has introduced a new particle species
	if(part)
	  AddV2(part, Find(part->GetPdgCode(), kTRUE)->GetV2(part->Pt()), fluc);
    }
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::GenerateOnTheFlyTracks(Int_t mult, Int_t pid, TClonesArray* event, Double_t fluc) 
{
    // generate tracks for the event. the tracks that are generated here are the mothers
    Double_t ph, et, pt, mm, px, py, pz, ee;
    // access the descired generator through the static find method. if the particle species is unknown, we add it to the generator
    NaiveFlowAndSpectrumGenerator* generator = Find(pid, kTRUE);
    Int_t _tempCounter = OnTheFlyEventGeneratorCounter;
    for(int i=_tempCounter; i<=mult+_tempCounter; ++i) {
        OnTheFlyEventGeneratorCounter++;
        TParticle* part = (TParticle*)event->ConstructedAt(i);
        part->SetStatusCode(1);
        part->SetFirstMother(-1);
        part->SetLastMother(-1);
        part->SetFirstDaughter(-1);
        part->SetLastDaughter(-1);
        ph = gRandom->Uniform(0,TMath::TwoPi());
        et = gRandom->Uniform(-1.0,+1.0);
        pt = generator->GetPt();
        part->SetPdgCode( pid );
        mm =  part->GetMass();
        px = pt*TMath::Cos(ph);
        py = pt*TMath::Sin(ph);
        pz = pt*TMath::SinH(et);
        ee = TMath::Sqrt( mm*mm + pt*pt+pz*pz );
        part->SetMomentum( px, py, pz, ee );
        // if requested add v2 to the sample of mothers
        if(fAddV2Mothers) AddV2(part, generator->GetV2(pt), fluc);
    }
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::DecayOnTheFlyTracks(TClonesArray *event) 
{
    // decay the tracks using a decayer that is set in the GenerateEventsOnTheFly.C macro
    if(!fDecayer) return;       // shouldn't happen ...
    Int_t kf;
    TClonesArray*  arr = new TClonesArray("TParticle",10);
    Int_t secondaries=0;
    Int_t nStart=0;
    Int_t nTracks = event->GetEntriesFast();
    TParticle *part, *part1;
    TClonesArray &in = *event;
    for(Int_t loop=0; loop!=3; ++loop) {
        secondaries=0;
        for(Int_t nTrack=nStart; nTrack!=nTracks; ++nTrack) {
            arr->Clear();
            part = (TParticle*) event->At( nTrack );
            if(!part) continue;
            kf = TMath::Abs(part->GetPdgCode());
            if(kf==22 )  ForceGammaDecay(arr, part);     // if gamma, we decay it by hand
            else {      // else we let pythia do it
                TLorentzVector pmom( part->Px(), part->Py(), part->Pz(), part->Energy() );
                fDecayer->Decay(kf,&pmom);
                fDecayer->ImportParticles(arr);
            }
            Int_t nDaughters = arr->GetEntries();
            if(nDaughters<2) continue;
            for(Int_t nDaughter=1; nDaughter!=nDaughters; ++nDaughter) {
                part1 = (TParticle*) (arr->At(nDaughter));
                if(!part1) continue; // was part1, mistake?
	        new(in[nTracks+secondaries]) TParticle( part1->GetPdgCode(),
						        part1->GetStatusCode(),
						        part1->GetFirstMother(),
						        part1->GetSecondMother(),
						        part1->GetFirstDaughter(),
						        part1->GetLastDaughter(),
						        part1->Px(),
						        part1->Py(),
						        part1->Pz(),
						        part1->Energy(),
						        part1->Vx(),
						        part1->Vy(),
						        part1->Vz(),
						        part1->T());
	        secondaries++;
                if(nDaughter==1) part->SetFirstDaughter(nTracks+secondaries);
                else if ((nDaughters-1)==nDaughter) part->SetLastDaughter(nTracks+secondaries);
                //else part->SetDaughter(nDaughter,nTracks+secondaries);
                }
            }
        nStart = nTracks;
        nTracks = event->GetEntries();
    }
    delete arr;
    return;
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::ForceGammaDecay(TClonesArray* arr, TParticle* part)
{
    // force gama decay
    const Float_t mass1(gRandom->Uniform(0.0000001, .1));   // FIXME randomly sample photon 'mass'
    const Float_t mass_d1(.000511);                         // daughter mass (electron)
    const Float_t mass_d2(.000511);
    Double_t p_dec1=sqrt((TMath::Abs(mass1*mass1-(mass_d1+mass_d2)*(mass_d1+mass_d2))*(mass1*mass1-(mass_d1-mass_d2)*(mass_d1-mass_d2))))/2/mass1;      // energy
    TLorentzVector *p_d1=new TLorentzVector();              // lorentz vectors for daughters
    TLorentzVector *p_d2=new TLorentzVector();
    Double_t yy_m=part->Y();                                // pseudo rapidity of mother
    Double_t pt_m=part->Pt();                               // pt of mother
    Double_t mt_m=sqrt(mass1*mass1+pt_m*pt_m);              // transverse mass of mother
    Double_t pz_m=mt_m*TMath::SinH(yy_m);                   // MAGIC
    Double_t phi_m=gRandom->Rndm()*2*TMath::Pi();
    Double_t px_m=pt_m*TMath::Cos(phi_m);
    Double_t py_m=pt_m*TMath::Sin(phi_m);
    Double_t costh_d=2*gRandom->Rndm()-1;
    Double_t phi_d=gRandom->Rndm()*2*TMath::Pi();
    Double_t pz_d=p_dec1*costh_d;
    Double_t pt_d=p_dec1*sqrt(1-costh_d*costh_d);
    Double_t px_d=pt_d*TMath::Cos(phi_d);
    Double_t py_d=pt_d*TMath::Sin(phi_d);
    p_d1->SetPxPyPzE(px_d,py_d,pz_d,sqrt(mass_d1*mass_d1+p_dec1*p_dec1));
    p_d2->SetPxPyPzE(-px_d,-py_d,-pz_d,sqrt(mass_d2*mass_d2+p_dec1*p_dec1));
    Double_t gamma_b=sqrt(mass1*mass1+pz_m*pz_m+pt_m*pt_m)/mass1;
    Double_t bx=px_m/gamma_b/mass1;    
    Double_t by=py_m/gamma_b/mass1;
    Double_t bz=pz_m/gamma_b/mass1;
    p_d1->Boost(bx,by,bz);
    p_d2->Boost(bx,by,bz);
    TParticle* daughter_1 = (TParticle*)arr->ConstructedAt(0);
    TParticle* daughter_2 = (TParticle*)arr->ConstructedAt(1);
    daughter_1->SetStatusCode(1);
    daughter_1->SetFirstMother(-1);
    daughter_1->SetLastMother(-1);
    daughter_1->SetFirstDaughter(-1);
    daughter_1->SetLastDaughter(-1);
    daughter_1->SetPdgCode(11);
    TLorentzVector& ref_p_d1 = *p_d1;
    daughter_1->SetMomentum(ref_p_d1);
    daughter_2->SetStatusCode(1);
    daughter_2->SetFirstMother(-1);
    daughter_2->SetLastMother(-1);
    daughter_2->SetFirstDaughter(-1);
    daughter_2->SetLastDaughter(-1);
    daughter_2->SetPdgCode(-11);
    TLorentzVector& pp = *p_d2;
    daughter_2->SetMomentum(pp);
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::InitGenerators()
{
   // initializes generators for a list of common particles
   // new generators will be created 'on the fly' if they are encountered
   // if a generator is requested for a non-existent pdg value, things will get unstable ...
   Short_t PDGcode[] = {11, -11, // e-               
                        12, // v_e
                        13, -13, // mu-
                        14, // v_mu
                        15, -15, // t-
                        16, // v_t
                        21, // g
                        22, // gamma
                        111, // pi0
                        211, -211, // pi+
                        113, // rho0
                        213, -213, // rho+
                        221, // eta
                        331, // eta prime
                        130, // k0l
                        310, // k0s
                        311, // k0
                        313, // k*
                        321, -321, // k+
                        323, -323, // k*+
                        223, // omega
                        411, -411, // d+
                        413, -413, // d*+
                        421, // d0
                        423, // d*0
                        333, // phi
                        443, // jpsi 
                        2112, // neutron
                        2212, -2212, // proton
                        3122, // lambda0
                        3312, -3312, // xi- (cascade)
                        3314, -3314, // xi*-
                        3322, // xi0
                        3324, // xi*0
                        3334}; // Omeg
   for(int i(0); i < 47; i++) fGenerators->Add(new NaiveFlowAndSpectrumGenerator(PDGcode[i], fQA, fFF));
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::PrintGenerators()
{
   // make sure all went ok
   if(!fGenerators) {
       printf(" > PANIC <\n");
       return;
   }
   printf(" > %i species available in generator\n", fGenerators->GetEntriesFast());
   for(int i(0); i < fGenerators->GetEntriesFast(); i++) 
       printf("   - PDG: %i \n", ((NaiveFlowAndSpectrumGenerator*)fGenerators->At(i))->GetPDGCode());
   printf(" more can be created 'on the fly', but beware of non-existent pdg values !");
}
//_____________________________________________________________________________
void AliFlowOnTheFlyEventGenerator::DoGeneratorQA(Bool_t v2, Bool_t v3) 
{
    // this function loops over all the species that are available in the fGenerators, which are
    // a) created by the user (in the event generator) as mother particles or
    // b) introduced by pythia (the decayer) as daughter particles
    // c) made by the InitGenerators() function
    // 'empty species' (i.e. species for which a generator exists but which were not actually sampled) are omitted
    // the output histograms are written to a file named "GeneratorfQA.root"
    if(!fQA) {
        printf(" > Request has been made for QA plots but QA histograms have been disabled !\n");
        return;
    }
    TFile *QAfile = new TFile("GeneratorfQA.root","RECREATE"); 
    for(int i(0); i < fGenerators->GetEntriesFast(); i++) {
        NaiveFlowAndSpectrumGenerator* gen = (NaiveFlowAndSpectrumGenerator*)fGenerators->At(i);
        if(!gen) continue;               // skip if we failed to get a generator
        TH1F* QApt = (TH1F*)gen->GetQAType(0); 
        TH2F* QAv2 = (TH2F*)gen->GetQAType(1);
        TH2F* QAv3 = (TH2F*)gen->GetQAType(2);
        if((!QApt)||(!QAv2)||(!QAv3)) {
            printf(" > couldn't read qa histogrmas for species %i <\n", gen->GetPDGCode());
            continue;
        }
        if(QApt->GetEntries()==0&&QAv2->GetEntries()==0&&QAv3->GetEntries()==0) continue; // skip if no tracks of this species have been sampled
        printf(" > saving QA histograms for sampled species %i <\n", gen->GetPDGCode());
        QAfile->mkdir(Form("fPDG_%i", gen->GetPDGCode()));        // create a directory in the output file
        QAfile->cd(Form("fPDG_%i", gen->GetPDGCode()));              // cd into this directory
        if(!QApt->GetEntries()==0) { // only plot the pt fSpectrum if this guy was generated
            if(!gen->GetPtSpectrum()->Integral(0,20)==0) {
                QApt->Scale(gen->GetPtSpectrum()->Integral(0,20)/QApt->Integral(),"width");
                QApt->Write();
            }
        }
        gen->GetPtSpectrum()->SetNpx(10000); // otherwise tf1 plot is very ugly
        gen->GetPtSpectrum()->Draw();
        gen->GetPtSpectrum()->Write();
        if(v2) {
            QAv2->Draw();
            gen->GetDifferentialV2()->SetNpx(10000);
            gen->GetDifferentialV2()->Draw();
            gen->GetDifferentialV2()->Write();
            QAv2->Write();
        }
        if(v3) {
            QAv3->Draw();
            gen->GetDifferentialV3()->Draw();
            gen->GetDifferentialV3()->SetNpx(10000);
            gen->GetDifferentialV3()->Write();
            QAv3->Write();
        }
    }
    QAfile->Close();
}
//_____________________________________________________________________________