ITS/UPGRADE/macros/PlotDoseExtrapolat.C

   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);
   6
   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");
  11 //
  12 //// Now the estimate of the occupancy
  13 //
  14 TF1 *occh = new TF1("Dose vs z at given r ",fdosez,-60.,60.,3);
  15 //
  16 Double_t r0=3.7;
  17 Double_t z0=14.1;
  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);
  23 tot/=2*z0;
  24 cout << tot << endl;
  25 // now you can normalize
  26 TF2 *occ;
  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);
  31 //occh->Draw();
  32 //c3->SetLogy();
  33 // 1 sigma
  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);
  41 occ->SetContour(20);
  42 occ->Draw("colz");
  43 occ->Write();
  44
  45 //
  46 f->Close();
  47 }
  48
  49
  50 Double_t DoseAtZAndR(Double_t z, Double_t r, Double_t occ0, Double_t z0, Double_t r0) {
  51 if(r<=0) {
  52   cout << "r should be positive " << endl;
  53  return 0;
  54 }
  55
  56 Double_t occ=0;
  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));
  60 cout << occ << endl;
  61 return TMath::Abs(occ);
  62 }
  63
  64
  65 Double_t focc(Double_t* zr, Double_t* par) {
  66 Double_t z=zr[0];
  67 Double_t r=zr[1];
  68 Double_t occ0=par[0];
  69 Double_t r0=par[1];
  70 return OccupancyAtZAndR(z,r,occ0,r0);
  71 }
  72
  73 Double_t fdosez(Double_t* zz, Double_t* par) {
  74 Double_t z=zz[0];
  75 Double_t r=par[2];
  76 Double_t pippo=par[0];
  77 Double_t r0=par[2];
  78 Double_t z0=par[1];
  79 return DoseAtZAndR(z,r,pippo,z0,r0);
  80 }
  81
  82
  83 Double_t fdose(Double_t* zr, Double_t* par) {
  84 Double_t z=zr[0];
  85 Double_t r=zr[1];
  86 Double_t dose0=par[0];
  87 Double_t z0=par[1];
  88 Double_t r0=par[2];
  89 return DoseAtZAndR(z,r,dose0,z0,r0);
  90 }
  91