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