TPC/Upgrade/macros/makeSpaceChargeMap.C

   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                        Double_t radialScaling = 2., Double_t epsilonScaling = 2./3.) {
   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 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;
  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
  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
  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 }