]>
Commit | Line | Data |
---|---|---|
60af5691 | 1 | // |
2 | // Origin: Christian Lippman, CERN, Christian.Lippmann@cern.ch | |
3 | // | |
4 | ||
5 | int makeSpaceChargeMap(Double_t multiplicity = 950., Double_t intRate = 5e4, Double_t eps = 10., | |
16206229 | 6 | Double_t gasfactor = 1., string filename = "SpaceChargeMap.root", |
7 | Double_t radialScaling = 2., Double_t epsilonScaling = 2./3.) { | |
60af5691 | 8 | // |
9 | // Charge distribution is splitted into two (RZ and RPHI) in order to speed up | |
10 | // the needed calculation time. It is dumped to | |
11 | // | |
12 | // Explanation of variables: | |
13 | // 1) multiplicity: charghed particle dn/deta for top 80% centrality (660 for 2011, | |
14 | // expect 950 for full energy) | |
15 | // 2) intRate: Total interaction rate (e.g. 50kHz for the upgrade) | |
16 | // 3) eps: Number of backdrifting ions per primary electron (0 for MWPC, e.g.10 for GEM) | |
17 | // 4) gasfactor: Use different gas. E.g. Ar/CO2 has twice the primary ionization, ion drift | |
18 | // velocity factor 2.5 slower, so gasfactor = 5. | |
19 | // | |
20 | ||
21 | TFile *f = new TFile(filename.c_str(), "RECREATE"); | |
22 | ||
23 | // some grid, not too coarse | |
24 | Int_t nr = 350; | |
25 | Int_t nphi = 180; | |
26 | Int_t nz = 500; | |
27 | ||
28 | const Double_t fgkIFCRadius= 83.5; // radius which renders the "18 rod manifold" | |
29 | const Double_t fgkOFCRadius= 254.5; // Mean Radius of the Outer Field Cage | |
30 | const Double_t fgke0 = 8.854187817e-12; // vacuum permittivity [A·s/(V·m)] | |
31 | ||
32 | Double_t dr = (fgkOFCRadius-fgkIFCRadius)/(nr+1); | |
33 | Double_t dphi = TMath::TwoPi()/(nphi+1); | |
34 | Double_t dz = 500./(nz+1); | |
35 | Double_t safty = 0.; // due to a root bug which does not interpolate the boundary .. | |
36 | // .. (first and last bin) correctly | |
37 | ||
38 | // Charge distribution in ZR (rotational symmetric) ------------------ | |
39 | ||
40 | TH2F *histoZR = new TH2F("chargeZR", "chargeZR", | |
41 | nr, fgkIFCRadius-dr-safty, fgkOFCRadius+dr+safty, | |
42 | nz, -250-dz-safty, 250+dz+safty); | |
43 | ||
44 | for (Int_t ir=1;ir<=nr;++ir) { | |
45 | Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir); | |
46 | for (Int_t iz=1;iz<=nz;++iz) { | |
47 | Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz); | |
48 | ||
49 | // recalculation to meter | |
50 | Double_t lZ = 2.5; // approx. TPC drift length | |
51 | Double_t rpM = rp/100.; // in [m] | |
52 | Double_t zpM = TMath::Abs(zp/100.); // in [m] | |
53 | ||
54 | // calculation of "scaled" parameters | |
55 | Double_t a = multiplicity*intRate/76628; | |
56 | //Double_t charge = gasfactor * ( a / (rpM*rpM) * (1 - zpM/lZ) ); // charge in [C/m^3/e0], no IBF | |
16206229 | 57 | Double_t charge = gasfactor * ( a / (TMath::Power(rpM,radialScaling)) * (1 - zpM/lZ + epsilonScaling*eps) ); // charge in [C/m^3/e0], with IBF |
60af5691 | 58 | |
59 | charge = charge*fgke0; // [C/m^3] | |
60 | ||
61 | // from MC simulation (Stefan) | |
62 | // for 50kHz | |
63 | Double_t kon = (2.62243e-09); // charge in [C/m^3] | |
64 | // Add to normal charge: gain 2000 with {0.25,0.5%) ion feedback | |
65 | //charge += eps*(kon/(rpM*rpM)); | |
66 | ||
67 | if (zp<0) charge *= 0.9; // Slightly less on C side due to front absorber | |
68 | ||
69 | histoZR->SetBinContent(ir, iz, charge); | |
70 | } | |
71 | } | |
72 | ||
73 | histoZR->Write("SpaceChargeInRZ"); | |
74 | ||
75 | // Charge distribution in RPhi (e.g. Floating GG wire) ------------ | |
76 | ||
77 | TH3F *histoRPhi = new TH3F("chargeRPhi", "chargeRPhi", | |
78 | nr, fgkIFCRadius-dr-safty, fgkOFCRadius+dr+safty, | |
79 | nphi, 0-dphi-safty, TMath::TwoPi()+dphi+safty, | |
80 | 2, -1, 1); // z part - to allow A and C side differences | |
81 | ||
82 | // some 'arbitrary' GG leaks | |
83 | Int_t nGGleaks = 5; | |
84 | Double_t secPosA[5] = {3,6,6,11,13}; // sector | |
85 | Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm | |
86 | Double_t secPosC[5] = {1,8,12,15,15}; // sector | |
87 | Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm | |
88 | ||
89 | for (Int_t ir=1;ir<=nr;++ir) { | |
90 | Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir); | |
91 | for (Int_t iphi=1;iphi<=nphi;++iphi) { | |
92 | Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi); | |
93 | for (Int_t iz=1;iz<=2;++iz) { | |
94 | Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz); | |
95 | ||
96 | Double_t charge = 0; | |
97 | ||
98 | for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks | |
99 | ||
100 | // A side | |
101 | Double_t secPos = secPosA[igg]; | |
102 | Double_t radialPos = radialPosA[igg]; | |
103 | ||
104 | if (zp<0) { // C side | |
105 | secPos = secPosC[igg]; | |
106 | radialPos = radialPosC[igg]; | |
107 | } | |
108 | ||
109 | // some 'arbitrary' GG leaks | |
110 | if ( (phip<(TMath::Pi()/9*(secPos+1)) && phip>(TMath::Pi()/9*secPos) ) ) { // sector slice | |
111 | if ( rp>(radialPos-2.5) && rp<(radialPos+2.5)) // 5 cm slice | |
112 | //charge = 300; | |
113 | charge = 0.; | |
114 | } | |
115 | ||
116 | } | |
117 | ||
118 | charge = charge*fgke0; // [C/m^3] | |
119 | histoRPhi->SetBinContent(ir,iphi,iz,charge); | |
120 | } | |
121 | } | |
122 | } | |
123 | ||
124 | histoRPhi->Write("SpaceChargeInRPhi"); | |
125 | ||
126 | f->Close(); | |
127 | ||
128 | } |