IonizationGar.cxx

Go to the documentation of this file.

#include "CgemDigitizerSvc/IonizationGar.h"
 
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/DataSvc.h"
 
#include "CLHEP/Units/PhysicalConstants.h"
 
#include <iomanip>
#include <iostream>
#include <fstream>
#include <cmath>
 
using namespace std;
using namespace Garfield;
 
IonizationGar::IonizationGar(){
}
 
IonizationGar::~IonizationGar(){
     delete m_gas;
     delete m_geo;
     delete m_comp;
     delete m_box;
     delete m_sensor;
}
 
void IonizationGar::init(unsigned int random, ICgemGeomSvc* geomSvc, double magConfig){
     m_random = random;
     m_geomSvc = geomSvc;
     m_magConfig = magConfig;
 
     initGeoGas();
}
 
void IonizationGar::setTrack(int particle, int charge, double p, double trkPosIn[], double trkPosOut[]){
     clear();
     int nE = 0;
 
     double x0 = 0., y0 = 0., z0 = 0., t0 = 0.;
 
     // //the length of track segment.
     // double d0 = sqrt(pow((trkPosOut[0]-trkPosIn[0]),2)+pow((trkPosOut[1]-trkPosIn[1]),2)+pow((trkPosOut[2]-trkPosIn[2]),2));
     // double d1 = sqrt(pow((trkPosOut[0]-trkPosIn[0]),2)+pow((trkPosOut[1]-trkPosIn[1]),2));
     // double trkL = d0;
 
     // double lambda = acos(d1/d0);
     // double theta = acos((trkPosOut[1]-trkPosIn[1])/d1);
       
     // double dx0 = cos(lambda)*sin(theta)/10.0; // mm->cm
     // double dy0 = cos(lambda)*cos(theta)/10.0; // mm->cm
     // double dz0 = sin(lambda)/10.0;         // mm->cm
 
     double dx0 = (trkPosOut[0] - trkPosIn[0])/10.0; // cm
     double dy0 = (trkPosOut[1] - trkPosIn[1])/10.0; // cm
     double dz0 = (trkPosOut[2] - trkPosIn[2])/10.0; // cm
     double trkL = sqrt(dx0*dx0 + dy0*dy0 + dz0*dz0);
 
     // cout << "track segment: " << setw(15) << trkPosIn[0] << setw(15) << trkPosIn[1] << setw(15) << trkPosIn[2]
     //       << setw(15) << trkPosOut[0] << setw(15) << trkPosOut[1] << setw(15) << trkPosOut[2] << endl;
     // cout << setw(15) << trkL << setw(15) << dx0 << setw(15) << dy0 << setw(15) << dz0 << endl;
 
     double minX = 0. < dx0 ? 0. : dx0;
     double maxX = 0. < dx0 ? dx0 : 0.;
     double minY = 0. < dy0 ? 0. : dy0;
     double maxY = 0. < dy0 ? dy0 : 0.;
     double minZ = 0. < dz0 ? 0. : dz0;
     double maxZ = 0. < dz0 ? dz0 : 0.;
     if(maxX<=minX){
      minX += -0.1;
      maxX += 0.1;
     }
     if(maxY<=minY){
      minY += -0.1;
      maxY += 0.1;
     }
     if(maxZ<=minZ){
      minZ += -0.1;
      maxZ += 0.1;
     }
 
     string particleType = getParticle(particle, charge);
 
     TrackHeed* track = new TrackHeed();
     track->SetSensor(m_sensor);
     if(m_debugging) track->EnableDebugging();
     else track->DisableDebugging();
     track->EnableMagneticField();
     track->SetParticle(particleType);
     track->SetMomentum(p*1.e9);
     // m_sensor->SetArea(-0.5, 0, -0.5, 0.5, 0.5, 0.5);
     m_sensor->SetArea(minX, minY, minZ, maxX, maxY, maxZ);
     // cout << setw(15) << minX << setw(15) << maxX << setw(15) << minY << setw(15) << maxY
     //       << setw(15) << minZ << setw(15) << maxZ << endl;
 
     // cout << "trkpos" << endl;
     // cout << setw(15) << trkPosIn[0] << setw(15) << trkPosIn[1] << setw(15) << trkPosIn[2]
     //       << setw(15) << trkPosOut[0] << setw(15) << trkPosOut[1] << setw(15) << trkPosOut[2] << endl;
 
     track->NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);
     // track->NewTrack(0, 0, 0, t0, 0, 1, 0);
     // cout << "particle: " << particleType << ", P = " << p << endl;
     // cout << setw(15) << x0 << setw(15) << y0 << setw(15) << z0 << setw(15) << t0
     //       << setw(15) << dx0 << setw(15) << dy0 << setw(15) << dz0 << endl;
 
     double xc = 0., yc = 0., zc = 0., tc = 0.; // cluster position
     int nc = 0;              // number of electrons
     double ec = 0.;              // energy
     double extra = 0.;               // not used
     int nCluster = 0;
 
     while(track->GetCluster(xc, yc, zc, tc, nc, ec, extra)){
          double length = sqrt(pow(xc,2)+pow(yc,2)+pow(zc,2));
      // cout << "length = " << length << ",  trkL = " << trkL << endl;
          if(length>trkL){
           cout << "overflow -------------------------" << endl;
           continue;
      }
      // cout << setw(5) << nCluster << setw(15) << xc << setw(15) << yc << setw(15) << zc << setw(15) << nc << endl;
          nCluster++;
          for(int j = 0; j<nc; ++j) {
               double xe, ye, ze, te;
               double ee, dxe, dye, dze;
               track->GetElectron(j, xe, ye, ze, te, ee, dxe, dye, dze); // unit: cm
 
               // cout << setw(5) << nE << setw(15) << xe << setw(15) << ye << setw(15) << ze << setw(15) << te << endl;
               m_ex.push_back(trkPosIn[0] + xe*10.0); // mm
               m_ey.push_back(trkPosIn[1] + ye*10.0); // mm
               m_ez.push_back(trkPosIn[2] + ze*10.0); // mm
               m_et.push_back(te);
               nE++;
          }
     }
     m_nIonE = nE;
     delete track;
}
 
string IonizationGar::getParticle(int particle, int charge) const{
     string particleName;
     if((0==particle)&&(-1==charge))      {particleName = "electron";} 
     else if((0==particle)&&(1==charge))  {particleName = "e+";}
     else if((1==particle)&&(1==charge))  {particleName = "mu+";}
     else if((1==particle)&&(-1==charge)) {particleName = "mu-";}
     else if((2==particle)&&(1==charge))  {particleName = "pi+";}
     else if((2==particle)&&(-1==charge)) {particleName = "pi-";}
     else if((3==particle)&&(1==charge))  {particleName = "K+";}
     else if((3==particle)&&(-1==charge)) {particleName = "K-";}
     else if((4==particle)&&(1==charge))  {particleName = "p";}
     else if((4==particle)&&(-1==charge)) {particleName = "p-bar";}
     else { 
      cout << "IonizationGar::getParticle() error"; 
      particleName = "error";
     }
     return particleName;
}
 
void IonizationGar::clear(){
     m_nIonE = 0;
     m_ex.clear();
     m_ey.clear();
     m_ez.clear();
     m_et.clear();
}
 
void IonizationGar::initGeoGas(){
     randomEngine.Seed(m_random);
     m_gas = new MediumMagboltz();
     m_gas->SetComposition("ar",90.,"isobutane",10.);
     m_gas->SetTemperature(293.15);
     m_gas->SetPressure(760.0);
     m_gas->SetMaxElectronEnergy(200.);
     m_gas->EnablePenningTransfer(0.44, 0.0, "Ar");
     m_gas->Initialise(true);
 
     m_box = new SolidBox(0., 0., 0., 10., 10., 10.); // center and half length
 
     m_geo = new GeometrySimple();
     m_geo->AddSolid(m_box, m_gas);
 
     m_comp = new ComponentAnalyticField();
     if(m_debugging) m_comp->EnableDebugging();
     else m_comp->DisableDebugging();
     m_comp->SetGeometry(m_geo);
 
     m_sensor = new Sensor();
     if(m_debugging) m_sensor->EnableDebugging();
     else m_sensor->DisableDebugging();
     m_sensor->AddComponent(m_comp);
}