Changeset 69 in Idarraga for allpix

Timestamp:

Feb 22, 2011, 5:41:33 PM (13 years ago)

Author:

benoit

Message:

Fixed Threshold value in FEI3 Digitizer

File:

• allpix/src/AllPixFEI3StandardDigitizer.cc (modified) (14 diffs)

allpix/src/AllPixFEI3StandardDigitizer.cc

@@ -28 +28 @@
 
 
-    //G4UnitDefinition ( name, symbol, category, value )
+        //G4UnitDefinition ( name, symbol, category, value )
 
         //// Unit for charge in FEI3
@@ -38 +38 @@
         AllPixGeoDsc * gD = GetDetectorGeoDscPtr();
 
-        m_digitIn.thl = 0;
+
+
+        m_digitIn.thl = 3500*elec;
 
         ///////////////////////////////////////////////////////
@@ -47 +49 @@
         //detectorThickness = 285*um;
 
-        electricField = biasVoltage/detectorThickness; // V/um
-        diffusionCoeff = 36.0*cm2/s; // [cm2/s]
-        chargeSharingConstant = 0.05;
+        electricFieldY = biasVoltage/detectorThickness; // V/um
+        electricFieldX = 0; // V/um
+        electricFieldZ = 0; // V/um
+
+
+        chargeSharingConstant = 0.02;
         fluence = 1e1/cm2;
         Beta_electrons = 5e-16*cm2/ns;
@@ -57 +62 @@
         // Geometry Related constants
         pitchX = gD->GetPixelX();
-        pitchY = gD->GetPixelY();
-        nPixX= gD->GetNPixelsX();
-        nPixY= gD->GetNPixelsY();
+        pitchY = gD->GetPixelY();
+        nPixX= gD->GetNPixelsX();
+        nPixY= gD->GetNPixelsY();
 
         //      pitchX = gD->GetSensorX();
         //      pitchY = gD->GetSensorY();
         //      nPixX= 18;
-        //nPixY= 164;
+        //nPixY= 164;
         ///////////////////////////////////////////////////////
+
+        // Silicon electron and hole transport constants
+
+        Temperature = 300.;
+
+
+        Default_Electron_Mobility=1400.0*cm2/s; // Electron mobility (cm2/Vs)
+        Default_Hole_Mobility=450.0*cm2/s;// Hole mobility (cm2/Vs
+
+        Default_Electron_D=36*cm2/s; // Electron mobility (cm2/s)
+        Default_Hole_D=12*cm2/s;// Hole mobility (cm2/s
+
+        //mobility dependence on electric field
+        Electron_AlphaField = 2.4e7*cm/s ; //[um/ns]
+        Electron_ThetaField = 0.8;
+        Electron_TempNominal = 600.0 ; // [K]
+
+        Electron_Beta = 2.0;
+        Electron_Saturation_Velocity = Electron_AlphaField*TMath::Power(1.0+Electron_ThetaField*TMath::Exp(Temperature/Electron_TempNominal),-1.0);
+
+        Hole_AlphaField = 2.4e7*cm/s ; //[um/ns]
+        Hole_ThetaField = 0.8;
+        Hole_TempNominal = 600.0 ; // [K]
+
+        Hole_Beta = 1.0 ;
+        Hole_Saturation_Velocity = Hole_AlphaField*TMath::Power(1.0+Hole_ThetaField*TMath::Exp(Temperature/Hole_TempNominal),-1.0);
+
+
         // physics switches
         doTrapping =false;
+
+        //Numerical integration accuray
+        Target = 1e-4 ;
+        tlow =  0.001*ns;
+        tup=    0.1*ns;
+        dtIni = 0.01*ns;
 
 
@@ -90 +129 @@
 }
 
+
+vector<G4double>  AllPixFEI3StandardDigitizer::ComputeDriftTimeFullField(G4double x, G4double y, G4double z)
+{
+        G4double driftTime=0;
+        G4double dt=dtIni;
+        G4double xtemp=x;
+        G4double ytemp=y;
+        G4double ztemp=z;
+
+        vector<G4double> step;
+
+        while(z>0 && z<detectorThickness){
+
+                driftTime+=dt;
+                step=RKF5Integration(x,y,z,dt);
+                xtemp+=step[0];
+                ytemp+=step[1];
+                ztemp+=step[2];
+                SetDt(dt,step[3]);
+        }
+
+        vector<G4double> output;
+        output[0]=xtemp;
+        output[1]=ytemp;
+        output[2]=ztemp;
+        output[3]=driftTime;
+        return output;
+
+}
+
+
+void AllPixFEI3StandardDigitizer::ComputeElectricField(G4double x, G4double y, G4double z){}
+
+
+G4double AllPixFEI3StandardDigitizer::GetElectricFieldNorm(G4double x, G4double y, G4double z){
+
+        return TMath::Sqrt(electricFieldX*electricFieldX +electricFieldY*electricFieldY +electricFieldZ*electricFieldZ );
+}
+
+G4double  AllPixFEI3StandardDigitizer::MobilityElectron(G4double x, G4double y, G4double z){
+
+        ComputeElectricField(x,y,z);
+        G4double parElectricField = GetElectricFieldNorm(x,y,z);
+        G4double mobilite = Default_Electron_Mobility*
+        TMath::Power((1.0/(1.+ TMath::Power((Default_Electron_Mobility*parElectricField)/
+                        Electron_Saturation_Velocity,Electron_Beta))),1.0/Electron_Beta);
+        return mobilite;
+        //cout << mobilite << " " << parElectricField << endl;
+
+}
+
+
+G4double AllPixFEI3StandardDigitizer::MobilityHole(G4double x, G4double y, G4double z){
+
+        ComputeElectricField(x,y,z);
+        G4double parElectricField = GetElectricFieldNorm(x,y,z);
+        G4double mobilite = Default_Hole_Mobility*
+        TMath::Power((1.0/(1.+ TMath::Power((Default_Hole_Mobility*parElectricField)/
+                        Hole_Saturation_Velocity,Hole_Beta))),1.0/Hole_Beta);
+        return mobilite;
+
+
+}
+
+
 G4double AllPixFEI3StandardDigitizer::ComputeDriftTimeUniformField(AllPixTrackerHit *a_hit)
 {
@@ -95 +199 @@
         //G4cout << "Hit position with respect to pixel (um) "<< (a_hit->GetPosWithRespectToPixel().z() + detectorThickness/2.0)/um  <<endl;
         //return (a_hit->GetPosWithRespectToPixel().z() + detectorThickness/2.0 ) /(mobility*electricField);
-        return ( detectorThickness) /(mobility*electricField);
+        return ( detectorThickness) /(mobility*electricFieldY);
 
 }
@@ -101 +205 @@
 G4double AllPixFEI3StandardDigitizer::ComputeDiffusionRMS(G4double tDrift)
 {
-        return TMath::Sqrt(2.*diffusionCoeff*tDrift);
+        return TMath::Sqrt(2.*Default_Electron_D*tDrift);
 
 }
@@ -155 +259 @@
                 tempPixel.second = (*hitsCollection)[itr]->GetPixelNbY();
 
+/*
                 //G4cout << "[Digitizer] New Hit" << endl;
                 // post position of the hit
@@ -160 +265 @@
                 //double postPosY = (*hitsCollection)[itr]->GetPostPixelNbY();
                 //Fix FANO Factor !!
-                G4double eHitTruth=(*hitsCollection)[itr]->GetEdep();
+                //G4double eHitTruth=(*hitsCollection)[itr]->GetEdep();
+*/
+
+
                 G4double eHit = elec*CLHEP::RandGauss::shoot((*hitsCollection)[itr]->GetEdep()/elec,TMath::Sqrt((*hitsCollection)[itr]->GetEdep()/elec)*0.118);
-                
-                //pixelsContent[tempPixel] += eHit;
-
-                ///
-
-                G4double driftTime = ComputeDriftTimeUniformField((*hitsCollection)[itr]);
-                G4double sigma = ComputeDiffusionRMS(driftTime);
-                //G4cout << "Sigma : " << sigma/um << G4endl;
-                //G4cout << "Tdrift : " << driftTime/s << G4endl;
-                //G4cout << TString::Format("Sigma = %f, Drift time = %f",sigma/um,driftTime/ns) << endl;
-                // Let's visit all neighbours
-                //int x = tempPixel.first;
-                //int y = tempPixel.second;
-
                 G4double xpos = (*hitsCollection)[itr]->GetPosWithRespectToPixel().x();
                 G4double ypos = (*hitsCollection)[itr]->GetPosWithRespectToPixel().y();
+                G4double zpos = (*hitsCollection)[itr]->GetPosWithRespectToPixel().z();
+                G4double sigma;
+                G4double driftTime;
+                if(doFullField==false){
+                         driftTime = ComputeDriftTimeUniformField((*hitsCollection)[itr]);
+                         sigma = ComputeDiffusionRMS(driftTime);
+                }
+                else{
+                        vector<G4double> data = ComputeDriftTimeFullField(xpos,ypos,zpos);
+                        driftTime = data[3];
+                        sigma = ComputeDiffusionRMS(driftTime);
+                        xpos=data[0];
+                        ypos=data[1];
+
+                }
 
                 // see if need more energy in neighbor
@@ -184 +293 @@
 
                 //loop over all hit pixel neighbors
-                G4cout << TString::Format("[MC Truth] Hit Energy=%f, x=%f, y=%f pixel(%i,%i)",eHitTruth/elec,xpos/um,ypos/um,tempPixel.first,tempPixel.second) << endl;
+                //G4cout << TString::Format("[MC Truth] Hit Energy=%f, x=%f, y=%f pixel(%i,%i)",eHitTruth/elec,xpos/um,ypos/um,tempPixel.first,tempPixel.second) << endl;
                 //if(fabs(xpos)>=pitchX/2.-10*sigma && fabs(ypos)>pitchY/2.0-10*sigma){
-                if(fabs(ypos) >=pitchY/2.-3*sigma){
-                if( (fabs(xpos) >= pitchX/2.-3*sigma || fabs(ypos) >=pitchY/2.-3*sigma) ){
+                if(fabs(ypos) >=pitchY/2.-2*sigma){
+                if( (fabs(xpos) >= pitchX/2.-2*sigma || fabs(ypos) >=pitchY/2.-2*sigma) ){
                   // G4cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*****!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" <<endl ;
                   for(int i=-1;i<=1;i++){
@@ -197 +306 @@
                                         {
                                         //We compute contribution of the hit to each pixels
-                                          double Etemp = IntegrateGaussian(xpos/nm,ypos/nm,sigma/nm,(-pitchX/2.0 + i*pitchX)/nm,(-pitchX/2.+(i+1)*pitchX)/nm,(-pitchY/2 + j*pitchY)/nm,(-pitchY/2 + (j+1)*pitchY)/nm, eHit );
+
+                                        double Etemp = IntegrateGaussian(xpos/nm,ypos/nm,sigma/nm,(-pitchX/2.0 + i*pitchX)/nm,(-pitchX/2.+(i+1)*pitchX)/nm,(-pitchY/2 + j*pitchY)/nm,(-pitchY/2 + (j+1)*pitchY)/nm, eHit );
                                         if(doTrapping==true) pixelsContent[extraPixel]+=ApplyTrapping(driftTime,Etemp);
                                         else pixelsContent[extraPixel] +=Etemp;
 
-                                        G4cout << TString::Format("[Digitizer] Pixel %i %i Energy=%f, Energy after Trapping=%f",extraPixel.first,extraPixel.second,pixelsContent[extraPixel]/elec,ApplyTrapping(driftTime,pixelsContent[extraPixel])/elec) << endl;
-                                        //cout << TString::Format("Pixel %i %i, Energy collected = %f",extraPixel.first,extraPixel.second,Etemp) << endl;
+                                        //G4cout << TString::Format("[Digitizer] Pixel %i %i Energy=%f, Energy after Trapping=%f",extraPixel.first,extraPixel.second,pixelsContent[extraPixel]/elec,ApplyTrapping(driftTime,pixelsContent[extraPixel])/elec) << endl;
+                                        cout << TString::Format("Pixel %i %i, Energy collected = %f",extraPixel.first,extraPixel.second,Etemp/keV) << endl;
 
                                         };
@@ -212 +322 @@
                   if(doTrapping==true)pixelsContent[extraPixel] +=ApplyTrapping(driftTime,pixelsContent[tempPixel]);
                   else pixelsContent[extraPixel] +=eHit;
-                  G4cout << TString::Format("[Digitizer] Pixel %i %i Energy=%f, Energy after Trapping=%f",extraPixel.first,extraPixel.second,pixelsContent[extraPixel]/elec,ApplyTrapping(driftTime,pixelsContent[extraPixel])/elec) << endl;
+                  //G4cout << TString::Format("[Digitizer] Pixel %i %i Energy=%f, Energy after Trapping=%f",extraPixel.first,extraPixel.second,pixelsContent[extraPixel]/elec,ApplyTrapping(driftTime,pixelsContent[extraPixel])/elec) << endl;
 
                 };
@@ -296 +406 @@
                 // If the charge in a given pixel is over the threshold
                 double threshold =CLHEP::RandGauss::shoot(m_digitIn.thl,chipNoise);
+                cout << TString::Format("Threshold= %f", threshold/keV) << endl;
+                //if((*pCItr).second > m_digitIn.thl)
                 if((*pCItr).second > threshold)
                 {
@@ -328 +440 @@
 
 }
+
+
+
+
+
+vector<G4double>  AllPixFEI3StandardDigitizer::RKF5Integration(G4double x, G4double y, G4double z,G4double dt)
+{
+// This function transport using Euler integration, for field (Ex,Ey,Ez),
+// considered constant over time dt. The movement equation are those
+// of charges in semi-conductors, sx= mu*E*dt;;
+        double k1x,k2x,k3x,k4x,k5x,k6x;
+        double k1y,k2y,k3y,k4y,k5y,k6y;
+        double k1z,k2z,k3z,k4z,k5z,k6z;
+        double dx,dy,dz;
+
+      ComputeElectricField(x,y,z);
+
+      k1x=-MobilityElectron(x,y,z)*electricFieldX*dt;
+      k1y=-MobilityElectron(x,y,z)*electricFieldY*dt;
+      k1z=-MobilityElectron(x,y,z)*electricFieldZ*dt;
+
+      ComputeElectricField(x+k1x/4,y+k1y/4,z+k1z/4);
+
+      k2x=-MobilityElectron(x+k1x/4,y+k1y/4,z+k1z/4)*electricFieldX*dt;
+      k2y=-MobilityElectron(x+k1x/4,y+k1y/4,z+k1z/4)*electricFieldY*dt;
+      k2z=-MobilityElectron(x+k1x/4,y+k1y/4,z+k1z/4)*electricFieldZ*dt;
+
+     ComputeElectricField(x+(9./32)*k2x+(3./32)*k1x,y+(9./32)*k2y+(3./32)*k1y,z+(9./32)*k2z+(3./32)*k1z);
+
+      k3x=-MobilityElectron(x+(9./32)*k2x+(3./32)*k1x,y+(9./32)*k2y+(3./32)*k1y,z+(9./32)*k2z+(3./32)*k1z)*electricFieldX*dt;
+      k3y=-MobilityElectron(x+(9./32)*k2x+(3./32)*k1x,y+(9./32)*k2y+(3./32)*k1y,z+(9./32)*k2z+(3./32)*k1z)*electricFieldY*dt;
+      k3z=-MobilityElectron(x+(9./32)*k2x+(3./32)*k1x,y+(9./32)*k2y+(3./32)*k1y,z+(9./32)*k2z+(3./32)*k1z)*electricFieldZ*dt;
+
+      ComputeElectricField(x-(7200./2197)*k2x+(1932./2197)*k1x+(7296./2197)*k3x,y-(7200./2197)*k2y+(1932./2197)*k1y+(7296./2197)*k3y,z-(7200./2197)*k2z+(1932./2197)*k1z+(7296./2197)*k3z);
+
+      k4x=-MobilityElectron(x-(7200./2197)*k2x+(1932./2197)*k1x+(7296./2197)*k3x,y-(7200./2197)*k2y+(1932./2197)*k1y+(7296./2197)*k3y,z-(7200./2197)*k2z+(1932./2197)*k1z+(7296./2197)*k3z)*electricFieldX*dt;
+      k4y=-MobilityElectron(x-(7200./2197)*k2x+(1932./2197)*k1x+(7296./2197)*k3x,y-(7200./2197)*k2y+(1932./2197)*k1y+(7296./2197)*k3y,z-(7200./2197)*k2z+(1932./2197)*k1z+(7296./2197)*k3z)*electricFieldY*dt;
+      k4z=-MobilityElectron(x-(7200./2197)*k2x+(1932./2197)*k1x+(7296./2197)*k3x,y-(7200./2197)*k2y+(1932./2197)*k1y+(7296./2197)*k3y,z-(7200./2197)*k2z+(1932./2197)*k1z+(7296./2197)*k3z)*electricFieldZ*dt;
+
+      ComputeElectricField(x-(8)*k2x+(439./216)*k1x+(3680./513)*k3x-(845./4104)*k4x,y-(8)*k2y+(439./216)*k1y+(3680./513)*k3y-(845./4104)*k4y,z-(8)*k2z+(439./216)*k1z+(3680./513)*k3z-(845./4104)*k4z);
+
+      k5x=-MobilityElectron(x-(8)*k2x+(439./216)*k1x+(3680./513)*k3x-(845./4104)*k4x,y-(8)*k2y+(439./216)*k1y+(3680./513)*k3y-(845./4104)*k4y,z-(8)*k2z+(439./216)*k1z+(3680./513)*k3z-(845./4104)*k4z)*electricFieldX*dt;
+      k5y=-MobilityElectron(x-(8)*k2x+(439./216)*k1x+(3680./513)*k3x-(845./4104)*k4x,y-(8)*k2y+(439./216)*k1y+(3680./513)*k3y-(845./4104)*k4y,z-(8)*k2z+(439./216)*k1z+(3680./513)*k3z-(845./4104)*k4z)*electricFieldY*dt;
+      k5z=-MobilityElectron(x-(8)*k2x+(439./216)*k1x+(3680./513)*k3x-(845./4104)*k4x,y-(8)*k2y+(439./216)*k1y+(3680./513)*k3y-(845./4104)*k4y,z-(8)*k2z+(439./216)*k1z+(3680./513)*k3z-(845./4104)*k4z)*electricFieldZ*dt;
+
+      ComputeElectricField(x+(2)*k2x-(8./27)*k1x-(3544./2565)*k3x-(1859./4104)*k4x-(11./40)*k5x,
+      y+(2)*k2y-(8./27)*k1y-(3544./2565)*k3y-(1859./4104)*k4y-(11./40)*k5y,
+      z+(2)*k2z-(8./27)*k1z-(3544./2565)*k3z-(1859./4104)*k4z-(11./40)*k5z);
+
+      k6x=-MobilityElectron(x+(2)*k2x-(8./27)*k1x-(3544./2565)*k3x-(1859./4104)*k4x-(11./40)*k5x,
+              y+(2)*k2y-(8./27)*k1y-(3544./2565)*k3y-(1859./4104)*k4y-(11./40)*k5y,
+              z+(2)*k2z-(8./27)*k1z-(3544./2565)*k3z-(1859./4104)*k4z-(11./40)*k5z)*electricFieldX*dt;
+      k6y=-MobilityElectron(x+(2)*k2x-(8./27)*k1x-(3544./2565)*k3x-(1859./4104)*k4x-(11./40)*k5x,
+              y+(2)*k2y-(8./27)*k1y-(3544./2565)*k3y-(1859./4104)*k4y-(11./40)*k5y,
+              z+(2)*k2z-(8./27)*k1z-(3544./2565)*k3z-(1859./4104)*k4z-(11./40)*k5z)*electricFieldY*dt;
+      k6z=-MobilityElectron(x+(2)*k2x-(8./27)*k1x-(3544./2565)*k3x-(1859./4104)*k4x-(11./40)*k5x,
+              y+(2)*k2y-(8./27)*k1y-(3544./2565)*k3y-(1859./4104)*k4y-(11./40)*k5y,
+              z+(2)*k2z-(8./27)*k1z-(3544./2565)*k3z-(1859./4104)*k4z-(11./40)*k5z)*electricFieldZ*dt;
+
+      dx=((16./135)*k1x+(6656./12825)*k3x+(28561./56430)*k4x-(9./50)*k5x+(2./55)*k6x);
+      dy=((16./135)*k1y+(6656./12825)*k3y+(28561./56430)*k4y-(9./50)*k5y+(2./55)*k6y);
+      dz=((16./135)*k1z+(6656./12825)*k3z+(28561./56430)*k4z-(9./50)*k5z+(2./55)*k6z);
+
+      double Ex,Ey,Ez,Erreur;
+      Ex=((1./360)*k1x-(128./4275)*k3x-(2197./75240)*k4x-(1./50)*k5x+(2./55)*k6x);
+      Ey=((1./360)*k1y-(128./4275)*k3y-(2197./75240)*k4y-(1./50)*k5y+(2./55)*k6y);
+      Ez=((1./360)*k1z-(128./4275)*k3z-(2197./75240)*k4z-(1./50)*k5z+(2./55)*k6z);
+      Erreur=sqrt(Ex*Ex+Ey*Ey+Ez*Ez);
+
+      vector<G4double> newpoint;
+      newpoint[0]=(dx);
+      newpoint[1]=(dy);
+      newpoint[2]=(dz);
+      newpoint[3]=(Erreur);
+
+     return newpoint;
+
+}
+
+
+
+
+//______________________________________________________________________________
+G4double  AllPixFEI3StandardDigitizer::SetDt(G4double Dt, G4double ErreurMoy)
+{
+
+        if(isnan(ErreurMoy)){Dt=tup;}
+        else if(ErreurMoy > Target){ Dt*=0.9;}
+        else if(ErreurMoy < Target){ Dt*=1.1;};
+
+
+        if(Dt<tlow) Dt=tlow;
+        if(Dt>tup)  Dt=tup;
+
+        return Dt;
+
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+