ALIROOT-5769: Better handling of the error (R. Preghenella)
[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;
130aabbc 7 Int_t zz[100], nn[100], nalpha, NFrag, ztot, ntot, ndeu;
e0b1fdb4 8 void spectator(Float_t, Int_t*, Int_t*);
9
130aabbc 10 TH2F *hsp = new TH2F("hsp","Spectator protons vs b",50,0.,20.,90,0.,90.);
e0b1fdb4 11 hsp -> SetXTitle("b (fm)");
12 hsp -> SetYTitle("Num. of spectator protons");
13
130aabbc 14 TH2F *hsn = new TH2F("hsn","Spectator neutrons vs b",50,0.,20.,130,0.,130.);
e0b1fdb4 15 hsn -> SetXTitle("b (fm)");
16 hsn -> SetYTitle("Num. of spectator neutrons");
17
130aabbc 18 TH2F *hFragp = new TH2F("hFragp","Free spectator protons vs b",50,0.,20.,60,0.,60.);
e0b1fdb4 19 hFragp -> SetXTitle("b (fm)");
20 hFragp -> SetYTitle("Num. of free spectator protons");
21
130aabbc 22 TH2F *hFragn = new TH2F("hFragn","Free spectator neutrons vs b",50,0.,20.,80,0.,80.);
e0b1fdb4 23 hFragn -> SetXTitle("b (fm)");
24 hFragn -> SetYTitle("Num. of free spectator neutrons");
25
130aabbc 26 TF1 *funb = new TF1("funb","x",0,20);
e0b1fdb4 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();
130aabbc 41 nalpha = gallio->GetNalpha();
42 NFrag = gallio->GetFragmentNum();
e0b1fdb4 43//
130aabbc 44 // Attach neutrons
45 ztot = gallio->GetZtot();
46 ntot = gallio->GetNtot();
25bd4833 47 gallio->AttachNeutrons();
130aabbc 48 //
49 nspectpfree = (nspectp-ztot-2*nalpha);
50 nspectnfree = (nspectn-ntot-2*nalpha);
51
52 // Removing deuterons
53 ndeu = (Int_t) (nspectnfree*gallio->DeuteronNumber());
54 nspectpfree -= ndeu;
55 nspectnfree -= ndeu;
56 //
57 hFragp -> Fill(b, nspectpfree);
58 hFragn -> Fill(b, nspectnfree);
e0b1fdb4 59//
60// Print
61 if(debug == 1){
130aabbc 62 printf("\n b = %f",b);
63 printf(" #spect p = %d, #spect n = %d\n",nspectp,nspectn);
64 printf(" #spect p free = %d, #spect n free = %d\n",nspectpfree,nspectnfree);
65 printf(" #fragments = %f ",NFrag);
66 /*for(Int_t i=0; i<NFrag; i++){
e0b1fdb4 67 printf("\n ZZ[%d] = %d, NN[%d] = %d\n",i,zz[i],i,nn[i]);
130aabbc 68 }*/
69 printf(" NAlpha = %d, Ztot = %d, Ntot = %d\n\n",nalpha,ztot,ntot);
e0b1fdb4 70 }
71 delete gallio;
72 } //Event loop
130aabbc 73
74 TProfile *profsp = hsp->ProfileX("profsp",-1,-1,"s");
75 TProfile *profsn = hsn->ProfileX("profsn",-1,-1,"s");
76 TProfile *profFragp = hFragp->ProfileX("profFragp",-1,-1,"s");
77 TProfile *profFragn = hFragn->ProfileX("profFragn",-1,-1,"s");
e0b1fdb4 78
130aabbc 79
80 //***********************************************************
81 // #### ROOT initialization
82 gROOT->Reset();
83 gStyle->SetCanvasColor(10);
84 gStyle->SetFrameFillColor(10);
85 gStyle->SetPalette(1);
86 gStyle->SetOptStat(0);
87 //
88 TCanvas *c1 = new TCanvas("c1","Fragmentation I",0,0,700,700);
89 c1->cd();
90 c1->SetFillColor(kAzure+7);
91 c1->Divide(2,2);
92 c1->cd(1);
93 hsp -> Draw("colz");
94 c1->cd(2);
95 hsn -> Draw("colz");
96 c1->cd(3);
97 hFragp -> Draw("colz");
98 c1->cd(4);
99 hFragn -> Draw("colz");
100 //
101 c1->Print("Fragment1.gif");
102 //
103 TCanvas *c2 = new TCanvas("c2","Fragmentation II",300,10,700,700);
104 c2->cd();
105 c2->SetFillColor(kAzure+10);
106 c2->Divide(2,2);
107 c2->cd(1);
108 profsp ->SetLineColor(kRed); profsp ->SetLineWidth(2);
109 profsp -> Draw("colz");
110 c2->cd(2);
111 profsn ->SetLineColor(kRed); profsn ->SetLineWidth(2);
112 profsn -> Draw("colz");
113 c2->cd(3);
114 profFragp -> Draw("colz");
115 profFragp ->SetLineColor(kAzure); profFragp ->SetLineWidth(2);
116 c2->cd(4);
117 profFragn -> Draw("colz");
118 profFragn ->SetLineColor(kAzure); profFragn ->SetLineWidth(2);
119 //
120 c1->Print("Fragment2.gif");
e0b1fdb4 121
122}
123
124
125void spectator(Float_t b, Int_t* NSpectp, Int_t* NSpectn)
126{
127 Float_t SppvsB[6] = {3.633,-1.518,1.360,-.4899e-1,-.2398e-2,.1066e-3};
128 Float_t SpnvsB[6] = {5.639,-1.685,1.803,-.3129e-1,-.6618e-2,.2352e-3};
129 Float_t Sigmap[4] = {.5668,-.2200e-1,.3657e-3,-.2201e-5};
130 Float_t Sigman[4] = {.4185,-.9798e-2,.1052e-3,-.4238e-6};
131
132 Float_t rnsp = SppvsB[0]+SppvsB[1]*b+SppvsB[2]*(b*b)+SppvsB[3]*(b*b*b)+
133 SppvsB[4]*(b*b*b*b)+SppvsB[5]*(b*b*b*b*b);
134 Float_t rnsn = SpnvsB[0]+SpnvsB[1]*b+SpnvsB[2]*(b*b)+SpnvsB[3]*(b*b*b)+
135 SpnvsB[4]*(b*b*b*b)+SpnvsB[5]*(b*b*b*b*b);
136 Float_t snsp = Sigmap[0]+Sigmap[1]*rnsp+Sigmap[2]*(rnsp*rnsp)+Sigmap[3]*
137 (rnsp*rnsp*rnsp);
138 Float_t snsn = Sigman[0]+Sigman[1]*rnsn+Sigman[2]*(rnsn*rnsn)+Sigman[3]*
139 (rnsn*rnsn*rnsn);
140
141 snsp = snsp*rnsp;
142 snsn = snsn*rnsn;
143
144 Float_t xgaup = gRandom->Gaus(0.0,1.0);
145 snsp = snsp*xgaup;
146 Float_t xgaun = gRandom->Gaus(0.0,1.0);
147 snsn = snsn*xgaun;
148 rnsp=rnsp+snsp;
149 rnsn=rnsn+snsn;
150
151 *NSpectp = Int_t(rnsp);
152 *NSpectn = Int_t(rnsn);
153
154}