#include "PReaction.h" #include "PFireball.h" #include "PFixedDecay.h" #include "PParticle.h" #include "PChannel.h" #include "PUtils.h" #include "TRandom.h" #include "TString.h" #include Int_t sim_dilep(TString strmother="rho0", TString strdalitz="", Bool_t dimuon=kFALSE, Float_t Eb=25.F, Float_t T=0.160F, Float_t blast=0.4F, Int_t model=15, Double_t sigma=0.8, TString outDir="./", TString outFile="pluto.rho0.epem", Int_t nevts=10000) { //gSystem->Load("/hera/hades/user/ekrebs/pluto/pluto_v5.42.2/libPluto.so"); // Dilepton decays with Pluto // mother_particle -> l+ l- // or Dalitz decays // mother_particle -> dalitz_particle l+ l- // // Requires Pluto version 5.42.2 or higher // // Input parameters: // strmother: Pluto Id of mother particle // strdalitz: Pluto Id of additional particle in dalitz decays // dimuon: kTRUE: make dimuon decay, otherwise dielectron decay // Eb: beam energy in GeV // T: temperature in GeV // blast: radial expansion velocity in velocity/c // model: // 1st bit: Breit-Wigner distribution // 2nd bit: do VDM (~1/M^3) // 3rd bit: Boltzmann distribution // 4th bit: cutoff at dielectron/dimuon mass instead of 2 pion mass // sigma: width of gaussian dNdy distribution // outDir: output directory // outFile: name of output file // nevts: number of events gRandom->SetSeed(0); PUtils::SetSeed((Int_t)(gRandom->Rndm()*64000)); PFireball *source = new PFireball((char*)strmother.Data(), Eb, T, 0., 1., blast, 1., 0., 0., 0.); // Width of rapidity distribution set to match NA49 results. source->Print(); PParticle *mother = new PParticle(strmother.Data()); PParticle *dalitz = NULL; if(strdalitz != "" && strdalitz != "no") { dalitz = new PParticle(strdalitz.Data()); } PParticle *lp = NULL; PParticle *lm = NULL; PParticle *dilep = NULL; TString slp, slm, sll; // Select decay to dimuon or dielectron. if(dimuon) { lp = new PParticle("mu+"); lm = new PParticle("mu-"); dilep = new PParticle("dimuon"); slp = "mu+"; slm = "mu-"; sll = "mumu"; } else { lp = new PParticle("e+"); lm = new PParticle("e-"); dilep = new PParticle("dilepton"); slp = "e+"; slm = "e-"; sll = "ee"; } PParticle *s1[] = { source, mother }; PChannel *c1 = new PChannel(s1, 1, 1); PChannel *c2 = NULL; PChannel *c3 = NULL; PReaction *r = NULL; TString outputFileName = Form("%s/%s",outDir.Data(),outFile.Data()); //Add user-defined model without VDM M**3 scaling TString modelId(strmother+"_novdm_"+slm+"_"+slp); TString modelDe(strmother+" decay without VDM"); PFixedDecay *pmodel = new PFixedDecay(modelId.Data(), modelDe.Data(),-1); pmodel->Add(strmother+", parent"); pmodel->Add(slp+", daughter"); pmodel->Add(slm+", daughter"); // BW + no VDM + Boltzmann if((model & 1) && !(model & 2) && (model & 4)) { makeDistributionManager()->Add(pmodel); //This enables the model makeDistributionManager()->Enable(modelId.Data()); makeDistributionManager()->LinkDB(); } // Enable dielectron/dimuon cutoff instead of 2-pi cutoff. if(model & 8) { TString cutoffId = strmother + "_ee_"; if(dimuon) { cutoffId.Append("mu+_mu-"); } else { cutoffId.Append("e-_e+"); } makeDistributionManager()->Enable(cutoffId.Data()); makeDistributionManager()->LinkDB(); } else { source->SetSigma(sigma); } //Print out what we have: //makeDistributionManager()->Print("decay_models"); // BW only if((model & 1) && !(model & 2) && !(model & 4)) { TString command("loop:["+strmother+"]->SetM({"+strmother+"_bw}->SampleTotalMass())"); Double_t m = makeStaticData()->GetParticleMass(strmother.Data()); Double_t w = makeStaticData()->GetParticleTotalWidth(strmother.Data()); TString bwCutoff(Form("if ( ([%s]->M() > %.3f + %.3f * 4. ) || ( [%s]->M() < %.3f - %.3f * 4.) )", strmother.Data(), m, w, strmother.Data(),m, w)); command += "; " + bwCutoff + " goto loop"; c1->Do((char*)command.Data()); } // Options for PReaction. See comments further down for explations. Int_t f0 = 0; Int_t f1 = 0; Int_t f2 = 0; Int_t f3 = 0; // Define decay channels. if(dalitz) { if(mother->Is("Delta+") && dimuon) { makeStaticData()->AddDecay(-1, "Delta+ -> p + dimuon", "D+", "p,dimuon", 1.e-7); } Int_t nChannels = 3; PParticle *motherdecay[] = {mother, dalitz, dilep}; c2 = new PChannel(motherdecay, 2, 1, 1); PParticle *dilepdecay[] = {dilep,lp,lm}; c3 = new PChannel(dilepdecay, 2, 1, 1); PChannel *cc[] = {c1, c2, c3}; /*** cc: Double pointer to PChannel array. outputFileName: Name without suffix. nChannels: Size of PChannel array. f0: Output options for the ROOT file (default 0): 0: Only the tracked (i.e. stable) particles are stored in the ROOT le. 1: All the particles are stored in the ROOT file, including the composite particles. f1: Decay-mode options (default 0): obsolete f2: Vertex-calculation options (default 0): 0: This option is off (no vertex calculation). 1: Production vertices are calculated for those particles that are written on le (depend- ing on the output option). The origin is considered to be the parent, or beam and target vertex. This is also the case for the products of the rst channel. For particles pro- duced in subsequent channels the production vertex is calculated by adding straight- line segments successively, each of length obtained as the product of the parent vector velocity times a lifetime randomly sampled from an exponential decay-time distribu- tion. An assumption of absence of any external magnetic elds is implicit. f3: ASCII output options (default 0): 0: No ASCII output. 1: ASCII output files, formatted for input to HGeant5 . Irrespective of the output option, ASCII les contain only tracked particles. Invoked from the PDecayManager (see below) class, a separate ASCII le is opened for each reaction channel processed. 2: Common ASCII output le for all reaction channels processed by a PDecayManager. ***/ r = new PReaction(cc,(Char_t*)outputFileName.Data(),nChannels,f0,f1,f2,f3); } else { Int_t nChannels = 2; PParticle *motherdecay[] = {mother, lp, lm}; c2 = new PChannel(motherdecay, 2, 1, 1); PChannel *cc[] = {c1, c2}; r = new PReaction(cc,(Char_t*)outputFileName.Data(),nChannels,f0,f1,f2,f3); } r->Print(); r->loop(nevts); return 0; }