Updated macro
[u/mrichter/AliRoot.git] / ZDC / Fragment.C
CommitLineData
e0b1fdb4 1void Fragment(Int_t nev=0, Int_t debug=0)
2{
3 gROOT->Reset();
e0b1fdb4 4
5 Float_t b;
6 Int_t nspectp, nspectn, nspectpfree, nspectnfree;
7 Int_t zz[100], nn[100], nalpha, ztot, ntot;
8 void spectator(Float_t, Int_t*, Int_t*);
9
10 TH2F *hsp = new TH2F("hsp","Spectator protons vs b",100,0.,20.,100,0.,150.);
11 hsp -> SetXTitle("b (fm)");
12 hsp -> SetYTitle("Num. of spectator protons");
13
14 TH2F *hsn = new TH2F("hsn","Spectator neutrons vs b",100,0.,20.,100,0.,150.);
15 hsn -> SetXTitle("b (fm)");
16 hsn -> SetYTitle("Num. of spectator neutrons");
17
18 TH2F *hFragp = new TH2F("hFragp","Free spectator protons vs b",100,0.,20.,100,0.,150.);
19 hFragp -> SetXTitle("b (fm)");
20 hFragp -> SetYTitle("Num. of free spectator protons");
21
22 TH2F *hFragn = new TH2F("hFragn","Free spectator neutrons vs b",100,0.,20.,100,0.,150.);
23 hFragn -> SetXTitle("b (fm)");
24 hFragn -> SetYTitle("Num. of free spectator neutrons");
25
26 TF1 *funb = new TF1("funb","x",0,15);
27 for(Int_t ievent=0; ievent<nev; ievent++){
28 if(ievent%1000==0){printf("Processing event %d\n",ievent);}
29 b= Float_t(funb->GetRandom());
30 spectator(b, &nspectp, &nspectn);
31 hsp -> Fill(b,nspectp);
32 hsn -> Fill(b,nspectn);
33 AliZDCFragment *gallio = new AliZDCFragment(b);
34 for(Int_t j=0; j<=99; j++){
35 zz[j] =0;
36 nn[j] =0;
37 }
38//
39// Generation of fragments
25bd4833 40 gallio->GenerateIMF();
e0b1fdb4 41 Int_t NFrag = gallio->GetFragmentNum();
42//
43// Attach neutrons
44 ztot=0;
45 ntot=0;
25bd4833 46 gallio->AttachNeutrons();
e0b1fdb4 47 nspectpfree = nspectp-ztot-2*nalpha;
48 nspectnfree = nspectn-ntot-2*nalpha;
49 hFragp -> Fill(b,nspectpfree);
50 hFragn -> Fill(b,nspectnfree);
51//
52// Print
53 if(debug == 1){
54 printf("\n b = %f\n",b);
55 printf("\n #spect p = %d, #spect n = %d\n",nspectp,nspectn);
56 printf("\n #spect p free = %d, #spect n free = %d\n",nspectpfree,nspectnfree);
57 printf("\n #fragments = %f\n",NFrag);
58 for(Int_t i=0; i<NFrag; i++){
59 printf("\n ZZ[%d] = %d, NN[%d] = %d\n",i,zz[i],i,nn[i]);
60 }
61 printf("\n NAlpha = %d, Ztot = %d, Ntot = %d\n\n",nalpha,ztot,ntot);
62 }
63 delete gallio;
64 } //Event loop
65
66 TCanvas *c1 = new TCanvas("c1","Fragmentation",0,10,580,700);
67 c1->cd();
68 c1->SetFillColor(29);
69 c1->Divide(2,2);
70 c1->cd(1);
71 hsp -> Draw("box");
72 c1->cd(2);
73 hsn -> Draw("box");
74 c1->cd(3);
75 hFragp -> Draw("box");
76 c1->cd(4);
77 hFragn -> Draw("box");
78
79}
80
81
82void spectator(Float_t b, Int_t* NSpectp, Int_t* NSpectn)
83{
84 Float_t SppvsB[6] = {3.633,-1.518,1.360,-.4899e-1,-.2398e-2,.1066e-3};
85 Float_t SpnvsB[6] = {5.639,-1.685,1.803,-.3129e-1,-.6618e-2,.2352e-3};
86 Float_t Sigmap[4] = {.5668,-.2200e-1,.3657e-3,-.2201e-5};
87 Float_t Sigman[4] = {.4185,-.9798e-2,.1052e-3,-.4238e-6};
88
89 Float_t rnsp = SppvsB[0]+SppvsB[1]*b+SppvsB[2]*(b*b)+SppvsB[3]*(b*b*b)+
90 SppvsB[4]*(b*b*b*b)+SppvsB[5]*(b*b*b*b*b);
91 Float_t rnsn = SpnvsB[0]+SpnvsB[1]*b+SpnvsB[2]*(b*b)+SpnvsB[3]*(b*b*b)+
92 SpnvsB[4]*(b*b*b*b)+SpnvsB[5]*(b*b*b*b*b);
93 Float_t snsp = Sigmap[0]+Sigmap[1]*rnsp+Sigmap[2]*(rnsp*rnsp)+Sigmap[3]*
94 (rnsp*rnsp*rnsp);
95 Float_t snsn = Sigman[0]+Sigman[1]*rnsn+Sigman[2]*(rnsn*rnsn)+Sigman[3]*
96 (rnsn*rnsn*rnsn);
97
98 snsp = snsp*rnsp;
99 snsn = snsn*rnsn;
100
101 Float_t xgaup = gRandom->Gaus(0.0,1.0);
102 snsp = snsp*xgaup;
103 Float_t xgaun = gRandom->Gaus(0.0,1.0);
104 snsn = snsn*xgaun;
105 rnsp=rnsp+snsp;
106 rnsn=rnsn+snsn;
107
108 *NSpectp = Int_t(rnsp);
109 *NSpectn = Int_t(rnsn);
110
111}