]>
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); | |
7f1a1b08 | 43 | |
44 | // For the normalization to same integral as radial exponent = 2 | |
45 | Double_t radialExponent = -2.; // reference = 2 | |
46 | Double_t radiusInner = histoZR->GetXaxis()->GetBinCenter(1) / 100.;//in [m] | |
47 | Double_t radiusOuter = histoZR->GetXaxis()->GetBinCenter(nr) / 100.;//in [m] | |
48 | Double_t integralRadialExponent2 = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1) | |
49 | - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1); | |
50 | ||
51 | radialExponent = -radialScaling; // user set | |
52 | Double_t integralRadialExponentUser = 0.; | |
53 | if(radialScaling > 1 + 0.000001 || radialScaling < 1 - 0.000001 ) // to avoid n = -1 | |
54 | integralRadialExponentUser = TMath::Power(radiusOuter,radialExponent+1) * 1./(radialExponent+1) | |
55 | - TMath::Power(radiusInner,radialExponent+1) * 1./(radialExponent+1); | |
56 | else | |
57 | integralRadialExponentUser = TMath::Log(radiusOuter) - TMath::Log(radiusInner); | |
58 | ||
59 | Double_t normRadialExponent = integralRadialExponent2 / integralRadialExponentUser; | |
60af5691 | 60 | |
61 | for (Int_t ir=1;ir<=nr;++ir) { | |
62 | Double_t rp = histoZR->GetXaxis()->GetBinCenter(ir); | |
63 | for (Int_t iz=1;iz<=nz;++iz) { | |
64 | Double_t zp = histoZR->GetYaxis()->GetBinCenter(iz); | |
65 | ||
66 | // recalculation to meter | |
67 | Double_t lZ = 2.5; // approx. TPC drift length | |
68 | Double_t rpM = rp/100.; // in [m] | |
69 | Double_t zpM = TMath::Abs(zp/100.); // in [m] | |
70 | ||
71 | // calculation of "scaled" parameters | |
72 | Double_t a = multiplicity*intRate/76628; | |
73 | //Double_t charge = gasfactor * ( a / (rpM*rpM) * (1 - zpM/lZ) ); // charge in [C/m^3/e0], no IBF | |
7f1a1b08 | 74 | Double_t charge = normRadialExponent * gasfactor * ( a / (TMath::Power(rpM,radialScaling)) * (1 - zpM/lZ + epsilonScaling*eps) ); // charge in [C/m^3/e0], with IBF |
75 | ||
60af5691 | 76 | charge = charge*fgke0; // [C/m^3] |
77 | ||
78 | // from MC simulation (Stefan) | |
79 | // for 50kHz | |
80 | Double_t kon = (2.62243e-09); // charge in [C/m^3] | |
81 | // Add to normal charge: gain 2000 with {0.25,0.5%) ion feedback | |
82 | //charge += eps*(kon/(rpM*rpM)); | |
83 | ||
84 | if (zp<0) charge *= 0.9; // Slightly less on C side due to front absorber | |
85 | ||
86 | histoZR->SetBinContent(ir, iz, charge); | |
87 | } | |
88 | } | |
89 | ||
90 | histoZR->Write("SpaceChargeInRZ"); | |
91 | ||
92 | // Charge distribution in RPhi (e.g. Floating GG wire) ------------ | |
93 | ||
94 | TH3F *histoRPhi = new TH3F("chargeRPhi", "chargeRPhi", | |
95 | nr, fgkIFCRadius-dr-safty, fgkOFCRadius+dr+safty, | |
96 | nphi, 0-dphi-safty, TMath::TwoPi()+dphi+safty, | |
97 | 2, -1, 1); // z part - to allow A and C side differences | |
98 | ||
99 | // some 'arbitrary' GG leaks | |
100 | Int_t nGGleaks = 5; | |
101 | Double_t secPosA[5] = {3,6,6,11,13}; // sector | |
102 | Double_t radialPosA[5] = {125,100,160,200,230}; // radius in cm | |
103 | Double_t secPosC[5] = {1,8,12,15,15}; // sector | |
104 | Double_t radialPosC[5] = {245,120,140,120,190}; // radius in cm | |
105 | ||
106 | for (Int_t ir=1;ir<=nr;++ir) { | |
107 | Double_t rp = histoRPhi->GetXaxis()->GetBinCenter(ir); | |
108 | for (Int_t iphi=1;iphi<=nphi;++iphi) { | |
109 | Double_t phip = histoRPhi->GetYaxis()->GetBinCenter(iphi); | |
110 | for (Int_t iz=1;iz<=2;++iz) { | |
111 | Double_t zp = histoRPhi->GetZaxis()->GetBinCenter(iz); | |
112 | ||
113 | Double_t charge = 0; | |
114 | ||
115 | for (Int_t igg = 0; igg<nGGleaks; igg++) { // loop over GG leaks | |
116 | ||
117 | // A side | |
118 | Double_t secPos = secPosA[igg]; | |
119 | Double_t radialPos = radialPosA[igg]; | |
120 | ||
121 | if (zp<0) { // C side | |
122 | secPos = secPosC[igg]; | |
123 | radialPos = radialPosC[igg]; | |
124 | } | |
125 | ||
126 | // some 'arbitrary' GG leaks | |
127 | if ( (phip<(TMath::Pi()/9*(secPos+1)) && phip>(TMath::Pi()/9*secPos) ) ) { // sector slice | |
128 | if ( rp>(radialPos-2.5) && rp<(radialPos+2.5)) // 5 cm slice | |
129 | //charge = 300; | |
130 | charge = 0.; | |
131 | } | |
132 | ||
133 | } | |
134 | ||
135 | charge = charge*fgke0; // [C/m^3] | |
136 | histoRPhi->SetBinContent(ir,iphi,iz,charge); | |
137 | } | |
138 | } | |
139 | } | |
140 | ||
141 | histoRPhi->Write("SpaceChargeInRPhi"); | |
142 | ||
143 | f->Close(); | |
144 | ||
145 | } |