pA_cocktail_35.C

Go to the documentation of this file.

//TITLE Thermal source model of p+nucleus 3.5 AGeV 
 
{
    // To save space, only leptons branches are generated and written to file.
    // allows for different seeds each time the macro is started
    
    PUtils::SetSeed(0);
    
    // makeDistributionManager()->Disable("helicity_angles");
    
    Bool_t freeze=kFALSE;
    
    Float_t Eb    = 3.5;     // beam energy in AGeV
    Float_t ctr   = 1.0;     // 
    //Float_t bmax  = 3.9;     // max. impact parameter (corresponds to 60% of all)
    Float_t bmax  = 0.0;         // Poisson sampling
    if (bmax>0.) ctr = 1.;   // use b sampling instead of Poisson sampling
    
    Float_t T1    = 0.080;   // temperature in GeV (for pion 2-component spectra)
    Float_t T2    = 0.080;   // temperature in GeV (assume this for thermalized source)
    Float_t frac  = 1.0;     // fraction of pion low-T component (from Jehad's fit to QMD)
    Float_t blast = 0.0;     // radial expansion velocity
    
    Float_t A2    = 1.0;     // polar distribution (from KaoS pion data)
    Float_t A4    = 0.0;
    Float_t v1    = 0.0;     // side flow
    Float_t v2    = 0.0;     // elliptic flow
    
    // Multiplicties for 3.5 AGeV P+nucleus min. bias events
    Float_t Mprot  = 2.03*ctr;       // proton
    Float_t Mneut  = 2.03*ctr;       // neutron
    Float_t Mdeut  = 0.0*ctr;      // deuteron (= 10% of proton)
    // <Apart> = ctr*(8.9 + 8.9 + 2*0.89) = 19.6 = A/2
    Float_t Mpi0   = 0.60*ctr;      // pi0 multiplicity in p+C 3.5 AGeV
    Float_t Mpip   = 0.60*ctr;      // pi+
    Float_t Mpim   = 0.60*ctr;      // pi-
    Float_t Meta   = 0.031*ctr;     // eta
    Float_t Momega = 0.011*ctr;     // omega
    Float_t Mrho0  = 0.011*ctr;     // rho0  
    Float_t Mphi   = 0.0005*ctr;    // phi
    Float_t MDelta = 2.*3./2.*Mpi0; // Delta0 + Delta+  
    
    Float_t enhance = 500.;           // enhancement factor for mesons
    //Float_t enhance = 5000.; 
    // enhance = 1.;           // enhancement factor
    Float_t enhancepi = 4.;            //enhancement factor for pi
    //Float_t enhancepi = 1;            //enhancement factor for pi
    
    Mpi0   *= enhancepi;
    //Mpip *= enhancepi;
    Meta   *= enhance;
    Momega *= enhance;
    Mrho0  *= enhance;
    Mphi   *= enhance;
    MDelta *= enhance;
 
    PFireball *source1 = new PFireball("pi0",Eb,T1,T2,frac,blast,A2,A4,v1,v2);
    PFireball *source2 = new PFireball("eta",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
    PFireball *source3 = new PFireball("w",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
    PFireball *source4 = new PFireball("w",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
    
    // source4->getFuncE()->Draw();
    // return;
 
    PFireball *source5  = new PFireball("rho0",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
    PFireball *source6  = new PFireball("phi",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
    PFireball *source7  = new PFireball("phi",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
    PFireball *source8  = new PFireball("D0",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
    PFireball *source9  = new PFireball("phi",Eb,T1,T2,frac,blast,A2,A4,v1,v2);
    PFireball *source10 = new PFireball("pi-",Eb,T1,T2,frac,blast,A2,A4,v1,v2);
    PFireball *source11 = new PFireball("p",Eb,T1,T2,frac,blast,A2,A4,v1,v2);
    PFireball *source12 = new PFireball("eta",Eb,T2,0.,1.,blast,0.0,A4,v1,v2);
 
    source1->SetW(1.0/enhancepi);
    source2->SetW(1.0/enhance);
    source3->SetW(1.0/enhance);
    source4->SetW(1.0/enhance);
    source5->SetW(1.0/enhance);
    source6->SetW(1.0/enhance);
    source7->SetW(1.0/enhance);
    source8->SetW(1.0/enhance);
    source9->SetW(1.0);
    source10->SetW(1.0);
    source11->SetW(1.0);
    source12->SetW(1.0/enhance);
    
    
    PChannel *c1 = source1->makeChannel(20, 0.012*Mpi0);          // pi0 decay directly
    PChannel *c2 = source2->makeChannel(1, 0.006*Meta);     // eta -> gamma e+e-
    PChannel *c3 = source3->makeChannel(1, 5.9e-4*Momega);  // omega -> pi0 e+e-
    PChannel *c4 = source4->makeChannel(1, 7.15e-5*Momega); // omega -> e+e-
    PChannel *c5 = source5->makeChannel(1, 4.48e-5*Mrho0);  // rho0 -> e+e-
    PChannel *c6 = source6->makeChannel(1, 3.00e-4*Mphi);   // phi -> e+e-
    PChannel *c7 = source7->makeChannel(1, 1.30e-4*Mphi);   // phi -> eta e+e-
    PChannel *c8 = source8->makeChannel(1, 4.4e-5*MDelta);  // Delta -> N e+e-
    PChannel *c9 = source9->makeChannel(1,0.49*Mphi); //phi->K+K-
    /*
      PChannel *c9 = source9->makeChannel(10, Mpip);  // Pion plus
      PChannel *c10 = source10->makeChannel(10, Mpim);  // Pion minus
      PChannel *c11 = source11->makeChannel(10, Mprot);  // Proton
    */
    PChannel *c12 = source12->makeChannel(1, 0.21*Meta);     // eta -> gamma, gamma 
    //PChannel *c14 = source13->makeChannel(1, 0.99*Mpi0);     // pi0 -> gamma, gamma
    //
    
    source1->Print();
    source2->Print();
    source3->Print();
    source4->Print();
    source5->Print();
    source6->Print();
    source7->Print();
    source8->Print();
    source9->Print();
    source10->Print();
    source11->Print();
    source12->Print();
    
    
    PParticle **ptcls;
    
    
    PParticle *gam1  = new PParticle("g");
    PParticle *diel1 = new PParticle("dilepton");
    PParticle *elec1 = new PParticle("e-");
    PParticle *posi1 = new PParticle("e+");
    
    
    ptcls = c1->GetParticles();
    PParticle *s1_1[] = {ptcls[1], gam1, diel1};
    PParticle *s1_2[] = {diel1, elec1, posi1};
    PChannel  *c1_1 = new PChannel(s1_1,2,1);
    PChannel  *c1_2 = new PChannel(s1_2,2,1);

    
    PParticle *gam2  = new PParticle("g");
    PParticle *diel2 = new PParticle("dilepton");
    PParticle *elec2 = new PParticle("e-");
    PParticle *posi2 = new PParticle("e+");

    
    ptcls = c2->GetParticles();
    PParticle *s2_1[] = {ptcls[1], gam2, diel2};
    PParticle *s2_2[] = {diel2, elec2, posi2};
    PChannel  *c2_1 = new PChannel(s2_1,2,1);
    PChannel  *c2_2 = new PChannel(s2_2,2,1);
    
    
    
    PParticle *pi03  = new PParticle("pi0");
    PParticle *diel3 = new PParticle("dilepton");
    PParticle *elec3 = new PParticle("e-");
    PParticle *posi3 = new PParticle("e+");
    
    ptcls = c3->GetParticles();
    PParticle *s3_1[] = {ptcls[1], pi03, diel3};
    PParticle *s3_2[] = {diel3, elec3, posi3};
    PChannel  *c3_1 = new PChannel(s3_1,2,1);
    PChannel  *c3_2 = new PChannel(s3_2,2,1);

    
    PParticle *elec4 = new PParticle("e-");
    PParticle *posi4 = new PParticle("e+");
 
    ptcls = c4->GetParticles();
    PParticle *s4_1[] = {ptcls[1], elec4, posi4};
    PChannel  *c4_1 = new PChannel(s4_1,2,1);

    
    PParticle *elec5 = new PParticle("e-");
    PParticle *posi5 = new PParticle("e+");
    
    ptcls = c5->GetParticles();
    PParticle *s5_1[] = {ptcls[1], elec5, posi5};
    PChannel  *c5_1 = new PChannel(s5_1,2,1);

    
    PParticle *elec6 = new PParticle("e-");
    PParticle *posi6 = new PParticle("e+");
    
    ptcls = c6->GetParticles();
    PParticle *s6_1[] = {ptcls[1], elec6, posi6};
    PChannel  *c6_1 = new PChannel(s6_1,2,1);

 
    PParticle *eta7  = new PParticle("eta");
    PParticle *diel7 = new PParticle("dilepton");
    PParticle *elec7 = new PParticle("e-");
    PParticle *posi7 = new PParticle("e+");
 
    ptcls = c7->GetParticles();
    PParticle *s7_1[] = {ptcls[1], eta7, diel7};
    PParticle *s7_2[] = {diel7, elec7, posi7};
    PChannel  *c7_1 = new PChannel(s7_1,2,1);
    PChannel  *c7_2 = new PChannel(s7_2,2,1);
    
    
    PParticle *neut8 = new PParticle("n");
    PParticle *diel8 = new PParticle("dilepton");
    PParticle *elec8 = new PParticle("e-");
    PParticle *posi8 = new PParticle("e+");
 
    ptcls = c8->GetParticles();
    PParticle *s8_1[] = {ptcls[1], neut8, diel8};
    PParticle *s8_2[] = {diel8, elec8, posi8};
    PChannel  *c8_1 = new PChannel(s8_1,2,1);
    PChannel  *c8_2 = new PChannel(s8_2,2,1);
 
    //phi->K+,k- decay
    
    ptcls= c9->GetParticles();
    PParticle *kp    = new PParticle("K+");
    PParticle *km    = new PParticle("K-");
    PParticle *s9_1[] = {ptcls[1],kp,km};
    PChannel  *c9_1 = new PChannel(s9_1,2,1);
    
    
    PParticle *gamma1 = new PParticle("g");
    PParticle *gamma2 = new PParticle("g");
 
    ptcls = c12->GetParticles();
    PParticle *s12_1[] = {ptcls[1], gamma1, gamma2};
    PChannel  *c12_1 = new PChannel(s12_1,2,1);

    
    
    
    
    PChannel *cc[] = {
    c1,c1_1,c1_2,        /* pi0          */
    c2,c2_1,c2_2,       /* eta Dalitz   */
    c3,c3_1,c3_2,       /* omega Dalitz */
    c4,c4_1,            /* omega direct */
    c5,c5_1,            /* rho0         */
    c6,c6_1,            /* phi direct   */
    c7,c7_1,c7_2,       /* phi Dalitz   */
    c8,c8_1,c8_2,        /* Delta Dalitz */
    //                 c9,c9_1,                  /* phi K+,k- production */
    //                  c10,                 /*Pim production*/
    //                  c11,                 /* Proton production */
    //                  c12,c12_1            /* eta 2 photon decay */
    };
    //change  *cc[]if you add want to add new channels and modify param after outFile
    
    PReaction *r = new PReaction(cc,"pA",21,1,0,0,1);
    //PReaction *r=new PReaction(cc,"pA",6,1,0,0,1);
    r->Print();
    
    //r->setUserSelection(selectLeptons);          // this works for a loaded function
    //r->setUserSelection(&select);                  //calls filter function in PReaction.cc
    //r->setTrigCond(1);    // require at least 1 lepton per event
    
    //r->SetDecayAll(1.);   // decay all particles with tau < 1 ns
    r->SetHGeant(0);      // set to 1, if PLUTO is run from HGeant prompt
    r->loop(10000);   // do n events
    
    //data->Draw("M()","(ID() == 51 || ID()==52 || ID()==41) * W()");
    //data->Draw("M()","(ID()==41) * W()","same"); 
    
}