1 Double_t OccupancyAtZAndR(Double_t z=0., Double_t r=2.2, Double_t occ0=2800, Double_t r0=3.7);
2 Double_t DoseAtZAndR(Double_t z=0., Double_t r=2.2, Double_t dose0=2800, Double_t z0=14.1, Double_t r0=3.7);
3 Double_t focc(Double_t* zr, Double_t* par);
4 Double_t fdosez(Double_t* zz, Double_t* par);
5 Double_t fdose(Double_t* zr, Double_t* par);
7 void PlotDoseExtrapolat(Bool_t zoom=kTRUE) {
8 //Double_t sigma=5.9; // PbPb 2010 run
9 Double_t sigma=7.94; // nominal parameter ( https://edms.cern.ch/file/445882/5/Vol_1_Chapter_21.pdf )
10 TFile *f = new TFile("prova.root","RECREATE");
12 //// Now the estimate of the occupancy
14 TF1 *occh = new TF1("Dose vs z at given r ",fdosez,-60.,60.,3);
18 Double_t dose0=2.2E+3 / 1.;
19 occh->SetParameter(0,1); // average occupancy at r0
20 occh->SetParameter(2,r0); // r0
21 occh->SetParameter(1,z0); // z0
22 Double_t tot=occh->Integral(-1.*z0,z0);
25 // now you can normalize
27 if (zoom) occ = new TF2("Dose vs z and r ",fdose,-20.,20.,1.5,5.,3); // zoomed
28 else occ = new TF2("Dose vs z and r ",fdose,-60.,60.,2.,50.,3); // not zoomed
29 TCanvas *c4=new TCanvas();
30 c4->SetRightMargin(0.15);
34 occ->SetParameter(0,dose0/tot); // occupancy at r0
35 occ->SetParameter(2,r0); // r0
36 occ->SetParameter(1,z0); // z0
37 occ->GetXaxis()->SetTitle("z (cm)");
38 occ->GetYaxis()->SetTitle("r (cm)");
39 occ->GetZaxis()->SetTitle("Dose (Gy)");
40 occ->GetZaxis()->SetTitleOffset(1.25);
50 Double_t DoseAtZAndR(Double_t z, Double_t r, Double_t occ0, Double_t z0, Double_t r0) {
52 cout << "r should be positive " << endl;
57 cout << " occ0=" << occ0 << " r0=" << r0 << " r=" << r << " z=" << z << endl;
58 if (TMath::Abs(z)<0.0000001) occ=occ0;
59 else occ=occ0*r0/r * TMath::Sin(TMath::ATan(r/z));
61 return TMath::Abs(occ);
65 Double_t focc(Double_t* zr, Double_t* par) {
70 return OccupancyAtZAndR(z,r,occ0,r0);
73 Double_t fdosez(Double_t* zz, Double_t* par) {
76 Double_t pippo=par[0];
79 return DoseAtZAndR(z,r,pippo,z0,r0);
83 Double_t fdose(Double_t* zr, Double_t* par) {
86 Double_t dose0=par[0];
89 return DoseAtZAndR(z,r,dose0,z0,r0);