TPC/SCDCorrections/sChargeMap.cxx

0001 #include <iostream>
0002 #include <fstream>
0003 #include "sChargeMap.h"
0004 #include "TMath.h"
0005 #include "TH2D.h"
0006 #include "TH3D.h"
0007 #include "TAxis.h"
0008 #include "TCanvas.h"
0009
0010 //=================
0011 sChargeMap::sChargeMap(int nr, float rmin, float rmax, int np, float pmin, float pmax, int nz, float hz, float ev, float iv) :
0012   fI(NULL),
0013   fE(NULL),
0014   fEcmperus(ev),
0015   fIcmperms(iv)
0016 {
0017   fI = new TH3D("fI","fI;r;phi;z",nr,rmin,rmax,np,pmin,pmax,nz,-hz,hz);
0018   fE = new TH3D("fE","fE;r;phi;z",nr,rmin,rmax,np,pmin,pmax,nz,-hz,hz);
0019 }
0020 //=================
0021 void sChargeMap::ScreenShot(char *name, char *ext, int n) {
0022   TCanvas *main = new TCanvas();
0023   main->Divide(2,2);
0024   TH2D *tmp = (TH2D*) fI->Project3D("xz");
0025   tmp->SetTitle(Form("IonMap (ev==%d)",n));
0026   main->cd(1)->SetLogz(1); tmp->Draw("colz");
0027   TH2D *tmp2 = (TH2D*) fE->Project3D("xz");
0028   tmp2->SetTitle(Form("ElectronMap (ev==%d)",n));
0029   main->cd(2)->SetLogz(1); tmp2->Draw("colz");
0030   TH1D *tmp1 = fI->ProjectionZ("IonMapZ");
0031   tmp1->SetTitle(Form("IonMapZ (ev==%d)",n));
0032   tmp1->Rebin(10);
0033   main->cd(3); tmp1->Draw("hist");
0034   TH1D *tmp12 = fE->ProjectionZ("ElectronMapZ");
0035   tmp12->SetTitle(Form("ElectronMapZ (ev==%d)",n));
0036   tmp12->Rebin(10);
0037   main->cd(4); tmp12->Draw("hist");
0038   main->SaveAs( Form("%s.%s",name,ext),ext);
0039   delete tmp;
0040   delete tmp1;
0041   delete tmp2;
0042   delete tmp12;
0043   delete main;
0044 }
0045 //=================
0046 void sChargeMap::SaveIonMap(char *name) {
0047   fI->SaveAs( name );
0048 }
0049 //=================
0050 void sChargeMap::SaveRho(char *name,int nrad,int nphi,int nhz) {
0051   int nr = fI->GetXaxis()->GetNbins();
0052   int np = fI->GetYaxis()->GetNbins();
0053   int nz = fI->GetZaxis()->GetNbins();
0054   nrad = TMath::Max( nr, nrad );
0055   nphi = TMath::Max( np, nphi );
0056   nhz  = TMath::Max( nz, nhz );
0057   TH3F *rho = new TH3F("rho","ChargeDensity [fC/cm^3];Radial [cm];Azimuthal [rad];Longitudinal [cm]",
0058                nrad, fI->GetXaxis()->GetBinLowEdge(1),fI->GetXaxis()->GetBinLowEdge( nr+1 ),
0059                nphi, 0,TMath::TwoPi(),
0060                nz,   fI->GetZaxis()->GetBinLowEdge(1),fI->GetZaxis()->GetBinLowEdge( nz+1 ));
0061   for(int i=0; i!=nrad; ++i) {
0062     float rmin = rho->GetXaxis()->GetBinLowEdge(i+1);
0063     float rmax = rho->GetXaxis()->GetBinLowEdge(i+2);
0064     int brmin = fI->GetXaxis()->FindBin( rmin+1e-6 );
0065     int brmax = fI->GetXaxis()->FindBin( rmax-1e-6 );
0066     for(int j=0; j!=nphi; ++j) {
0067       for(int k=0; k!=nhz; ++k) {
0068     float zmin = rho->GetZaxis()->GetBinLowEdge(k+1);
0069     float zmax = rho->GetZaxis()->GetBinLowEdge(k+2);
0070     int bzmin = fI->GetZaxis()->FindBin( zmin+1e-6 );
0071     int bzmax = fI->GetZaxis()->FindBin( zmax-1e-6 );
0072     float n = 0;
0073     for(int ii=brmin; ii!=brmax+1; ++ii)
0074       for(int jj=1; jj!=np+1; ++jj)
0075         for(int kk=bzmin; kk!=bzmax+1; ++kk)
0076           n += fI->GetBinContent( ii, jj, kk );
0077     float dv = TMath::Power(0.5*(rmax+rmin),2)*(rmax-rmin)*(zmax-zmin)*TMath::TwoPi()*(1./nphi);
0078     float drho = n / (dv*10000) * 1.602; // [fC/cm^3]
0079     //std::cout << "   n " << n << " || dv " << dv << std::endl;
0080     rho->SetBinContent( i+1, j+1, k+1, drho );
0081       }
0082     }
0083   }
0084   rho->SaveAs( name );
0085 }
0086 //=================
0087 void sChargeMap::Fill(float r, float p, float z, float w) {
0088   int vd1 = fI->GetZaxis()->FindBin(double(0));
0089   int vd2 = fI->GetZaxis()->FindBin(z);
0090   if(vd1==vd2) return;
0091   if(w<0) fE->Fill(r,p,z,TMath::Abs(w));
0092   else fI->Fill(r,p,z,w);
0093 }
0094 //=================
0095 void sChargeMap::Propagate(float ms) {
0096   int NZ = fI->GetZaxis()->GetNbins();
0097   //std::cout << " Ion Deltacm  " << fIcmperms*ms << std::endl;
0098   //std::cout << " Ion bin width  " << fI->GetZaxis()->GetBinWidth(1) << std::endl;
0099   int fIDis = fIcmperms*ms / fI->GetZaxis()->GetBinWidth(1);
0100   //std::cout << " DeltaBin  " << fIDis << std::endl;
0101   if(fIDis<1) {
0102     printf("Improve grid!!!!\n");
0103     fIDis=1;
0104   }
0105   // IONS: the easy part. Ions propagate to Anode plate (z==0) and disappear.
0106   // what i will do is to displace the bincontents according to the ion velocity
0107   // whenever the ions pass the Anode plate they will go to either under or overflow.
0108   int vd = fI->GetZaxis()->FindBin(double(0)); // anode
0109   // left-side
0110   for(int nz=vd-1; nz!=0; --nz) {
0111     int newz = nz+fIDis;
0112     if(newz>=vd) newz = 0; // to underflow
0113     for(int nr=0; nr!=fI->GetXaxis()->GetNbins(); ++nr)
0114       for(int np=0; np!=fI->GetYaxis()->GetNbins(); ++np) {
0115     float nc = fI->GetBinContent(nr+1,np+1,newz) + fI->GetBinContent(nr+1,np+1,nz);
0116     fI->SetBinContent(nr+1,np+1,nz, 0.0);
0117     fI->SetBinContent(nr+1,np+1,newz, nc );
0118       }
0119   }
0120   // right-side
0121   for(int nz=vd+1; nz!=NZ+1; ++nz) {
0122     int newz = nz-fIDis;
0123     if(newz<=vd) newz = NZ + 1; // to overflow
0124     for(int nr=0; nr!=fI->GetXaxis()->GetNbins(); ++nr)
0125       for(int np=0; np!=fI->GetYaxis()->GetNbins(); ++np) {
0126     float nc = fI->GetBinContent(nr+1,np+1,newz) + fI->GetBinContent(nr+1,np+1,nz);
0127     fI->SetBinContent(nr+1,np+1,nz,0.0);
0128     fI->SetBinContent(nr+1,np+1,newz, nc );
0129       }
0130   }
0131   //float fBF = 0.01; // percentage
0132   float fBF = 1; // percentage
0133   // ELECTRONS: though. Electrons will propagate according to the End plates (|z|==hz) and produce ions (backflow).
0134   // what i will do is to compute if the time passed is enough to reach the endplates in one go
0135   // if so, then the remaining is used to compute the position of ions backflowing and inject ions accordingly
0136   // if not, then the electrons are displaced normally
0137
0138   // left-side
0139   for(int nz=1; nz!=vd; ++nz) {
0140     float dz = fE->GetZaxis()->GetBinCenter(nz) - fE->GetZaxis()->GetBinLowEdge(1);
0141     float t2ep = dz / fEcmperus;
0142     float msres = ms - t2ep/1000;
0143     //std::cout << " Electron msres  " << msres << std::endl;
0144     if(msres>0) {
0145       // make backflow ions appear
0146       int newz = 1+fIcmperms*msres / fI->GetZaxis()->GetBinWidth(1);
0147       //std::cout << " Ion where in bin  " << newz << std::endl;
0148       if(newz>=vd) newz = 0; // to overflow
0149       for(int nr=0; nr!=fI->GetXaxis()->GetNbins(); ++nr)
0150     for(int np=0; np!=fI->GetYaxis()->GetNbins(); ++np) {
0151       float nc = fE->GetBinContent(nr+1,np+1,nz);
0152       if(nc>0) {
0153         nc = fI->GetBinContent(nr+1,np+1,newz) + fBF*fE->GetBinContent(nr+1,np+1,nz);
0154         //std::cout << " BACKFLOW: from " << nz << " to " << newz << " w=" << fE->GetBinContent(nr+1,np+1,nz) << " ==> " << nc << std::endl;
0155         fE->SetBinContent(nr+1,np+1,nz,0);
0156         fI->SetBinContent(nr+1,np+1,newz, nc );
0157       }
0158     }
0159     } else {
0160       // move electrons
0161       std::cout << " Electron Deltacm  " << fEcmperus*ms*1000 << std::endl;
0162       std::cout << " Electron bin width  " << fE->GetZaxis()->GetBinWidth(1) << std::endl;
0163       float fEDis = fEcmperus*ms*1000 / fE->GetZaxis()->GetBinWidth(1);
0164       std::cout << " DeltaBin  " << fEDis << std::endl;
0165       int newz = nz-fEDis;
0166       if(newz<1) newz = 0; // to underflow
0167       for(int nr=0; nr!=fI->GetXaxis()->GetNbins(); ++nr)
0168     for(int np=0; np!=fI->GetYaxis()->GetNbins(); ++np) {
0169       float nc = fE->GetBinContent(nr+1,np+1,newz) + fE->GetBinContent(nr+1,np+1,nz);
0170       fE->SetBinContent(nr+1,np+1,nz,0);
0171       fE->SetBinContent(nr+1,np+1,newz, nc );
0172     }
0173     }
0174   }
0175   // right-side
0176   for(int nz=vd+1; nz!=NZ+1; ++nz) {
0177     float dz = fE->GetZaxis()->GetBinLowEdge(NZ+1) - fE->GetZaxis()->GetBinCenter(nz);
0178     float t2ep = dz / fEcmperus;
0179     float msres = ms - t2ep/1000;
0180     //std::cout << " Electron msres  " << msres << std::endl;
0181     if(msres>0) {
0182       // make backflow ions appear
0183       int newz = NZ - fIcmperms*msres / fI->GetZaxis()->GetBinWidth(1);
0184       //std::cout << " Ion where in bin  " << newz << std::endl;
0185       if(newz<=vd) newz = NZ+1; // to overflow
0186       for(int nr=0; nr!=fI->GetXaxis()->GetNbins(); ++nr)
0187     for(int np=0; np!=fI->GetYaxis()->GetNbins(); ++np) {
0188       float nc = fE->GetBinContent(nr+1,np+1,nz);
0189       if(nc>0) {
0190         nc = fI->GetBinContent(nr+1,np+1,newz) + fBF*fE->GetBinContent(nr+1,np+1,nz);
0191         //std::cout << " BACKFLOW: from " << nz << " to " << newz << " w=" << fE->GetBinContent(nr+1,np+1,nz) << " ==> " << nc << std::endl;
0192         fE->SetBinContent(nr+1,np+1,nz,0);
0193         //fI->SetBinContent(nr+1,np+1,newz, nc );
0194       }
0195     }
0196     } else {
0197       // move electrons
0198       std::cout << " Electron Deltacm  " << fEcmperus*ms*1000 << std::endl;
0199       std::cout << " Electron bin width  " << fE->GetZaxis()->GetBinWidth(1) << std::endl;
0200       float fEDis = fEcmperus*ms*1000 / fE->GetZaxis()->GetBinWidth(1);
0201       std::cout << " DeltaBin  " << fEDis << std::endl;
0202       int newz = nz+fEDis;
0203       if(newz>NZ) newz = NZ+1; // to underflow
0204       for(int nr=0; nr!=fI->GetXaxis()->GetNbins(); ++nr)
0205     for(int np=0; np!=fI->GetYaxis()->GetNbins(); ++np) {
0206       float nc = fE->GetBinContent(nr+1,np+1,newz) + fE->GetBinContent(nr+1,np+1,nz);
0207       fE->SetBinContent(nr+1,np+1,nz,0);
0208       fE->SetBinContent(nr+1,np+1,newz, nc );
0209     }
0210     }
0211   }
0212 }