GARFDRIFT.h 5.33 KB
Newer Older
1 2 3 4 5 6 7 8 9 10
#include "Garfield/GeometrySimple.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewMedium.hh"
		 
11
double GARFDRIFT(double energy, double t_start, G4ThreeVector start, G4ThreeVector delta_pos, double R, int Pad_start, double ID,
12
		 double ScoLength, double SegLength, double event, double raw_energy_check) {
13

14 15 16 17 18 19
  ofstream file;
  /*        
  file.open("garf_E_check.dat", std::ios::app);
  file << raw_energy_check << "\t" << R << endl;
  file.close();
  */
20 21 22
  const double rWire = 1.2;
  const double rTube = 5.;
  const double lTube = 25.19;
23
  //TH1F *hist = new TH1F("","",1000,0,10000); // 10 ns bins
24 25 26 27 28
  //For R6 = 2, R7 = 1
  //const double vWire = -375.; //For V_total = 750V
  //const double vWire = -645.3312; //For V_total = 1200V //SHOULD BE 583.741 !?
  //const double vWire = -535.0959; //For V_total = 1100V
  //const double vWire = -875.6116; //For V_total = 1800V 
29 30 31
  //const double vWire =  -681.0312; //For V_total = 1400V
  //const double vWire = -730.; 
  const double vWire = -1000.;
32
  //const double vWire = -584.7139; // ND run 4230 (2008)
33 34
  const double vTube = 0.;
  double x_start, y_start, z_start;
35
  //double production_bias = 0.01;
36
  double production_bias = 0.001;
37
  double t_end, x_end, y_end, z_end;
38
  int status;
39 40
  double Pad_end;
  //double w_value = 35.; //for HeCO2 --guess
41 42
  //double w_value = 26.31; //for argon
  double w_value = 41.1;// for He
43 44 45
  energy = energy*production_bias;
  //int electrons  = (int)(energy/w_value*production_bias);
  int electrons = (int)(energy/w_value);
46
  double energy_excess = 0.;
47
  int sector;
48
  vector<double> pad_vec, time_vec; 
49 50 51 52 53 54 55 56 57
  
  Garfield::MediumMagboltz* gas = new Garfield::MediumMagboltz();
  Garfield::ViewMedium* mediumView = new Garfield::ViewMedium();
  Garfield::GeometrySimple* geo = new Garfield::GeometrySimple();
  Garfield::SolidTube* tube = new Garfield::SolidTube(0, 0, 0, 0, rTube, lTube/2);  
  Garfield::ComponentAnalyticField* cmp = new Garfield::ComponentAnalyticField();
  Garfield::Sensor* sensor = new Garfield::Sensor();
  Garfield::AvalancheMC* drift = new Garfield::AvalancheMC();
  
58 59 60 61 62 63 64 65 66
  //cout << "\n" << electrons << "\n" << endl;
  
  energy_excess = (energy - electrons*w_value)/production_bias;
  
  //  if(R>rWire) {

    //    energy_excess = (energy/w_value*production_bias - electrons)*w_value/production_bias;

    //    energy_excess = (energy/w_value - electrons)*w_value;
67 68
  
  if(R>rWire && electrons>0) {
69
    
70 71 72 73 74
      //gas->LoadGasFile("he_90_co2_10_200mbar.gas");
      //gas->LoadGasFile("ar_90_ch4_10.gas"); //atomic fraction in this case (P10).
      //gas->LoadGasFile("he_90_co2_10_500mbar.gas"); //mass fraction in this case
      //gas->LoadGasFile("he_90_co2_10_1bar.gas");

75
    gas->LoadGasFile("he_90_co2_10_100mbar.gas");
76
    //gas->LoadGasFile("he_90_co2_10_600mbar.gas");
77
    
78
    gas->Initialise(true);
79
      
80 81 82 83 84 85 86
    mediumView->SetMedium(gas);
    geo->AddSolid(tube, gas);
    cmp->SetGeometry(geo);
    cmp->AddWire(0, 0, 2*rWire, vWire, "c");
    cmp->AddTube(rTube, vTube, 0, "a");
    sensor->AddComponent(cmp);
    drift->SetSensor(sensor);
87 88 89

      //drift->DisableAttachment();
    
90 91 92
    for(int e=0;e<electrons;e++) {
      
      double randomize = (double)std::rand() / (double)RAND_MAX;
93
	
94 95 96
      x_start = start.x() + delta_pos.x()*randomize;
      y_start = start.y() + delta_pos.y()*randomize;
      z_start = start.z() + delta_pos.z()*randomize;
97 98 99 100 101
      /*
	file.open("garf_electrons_start_RZ.dat", std::ios::app);
	file <<  pow(pow(x_start,2)+pow(y_start,2),0.5) << "\t" << z_start <<  endl;	  
	file.close();
      */
102
      drift->DriftElectron(x_start, y_start, z_start, t_start);
103
      
104
      int np = drift->GetNumberOfDriftLinePoints();
105
      //drift->GetEndPoint(x_end, y_end, z_end, t_end, status);
106
      drift->GetDriftLinePoint(np-1, x_end, y_end, z_end, t_end); //np-1 is last DriftLineEndPoint
107 108
      //int np = drift->GetNuberOfElectronEndPoints();
      //drift->GetElectronEndPoint(np-1, x_end, y_end, z_end, t_end); 
109 110
      
      Pad_end = (int)((z_end + ScoLength / 2.) / SegLength ) + 1; //new Pad number 
111 112 113 114
      
      if(x_end > 0 and y_end > 0) {
	if(abs(x_end) > abs(y_end)) sector = 0;
	if(abs(x_end) < abs(y_end)) sector = 1;
115
      }
116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
      if(x_end < 0 and y_end > 0) {
	if(abs(x_end) < abs(y_end)) sector = 2;
	if(abs(x_end) > abs(y_end)) sector = 3;
      }
      if(x_end < 0 and y_end < 0) {
	if(abs(x_end) > abs(y_end)) sector = 4;
	if(abs(x_end) < abs(y_end)) sector = 5;
      }
      if(x_end > 0 and y_end < 0) {
	if(abs(x_end) < abs(y_end)) sector = 6;
	if(abs(x_end) > abs(y_end)) sector = 7;
      }
      
      
      //if(pow(pow(x_end,2)+pow(y_end,2),0.5)>4.9) { //reached the anode
      file.open("signal.dat", std::ios::app);
      file << event << "\t" << ID << "\t" << sector << "\t" <<  Pad_end << "\t" << t_end << "\t" << pow(pow(x_start,2)+pow(y_start,2),0.5) << "\t" << pow(pow(x_end,2)+pow(y_end,2),0.5) << endl;	  
      //file << event << "\t" << ID << "\t" << sector << "\t" <<  Pad_end << "\t" << t_end
      file.close();
      //}
      
137 138
    }
  }
139 140 141
  //}

delete gas; delete mediumView; delete geo; delete tube; delete cmp; delete sensor; delete drift; //otherwise these are held in memory and cause a killed: 9 crash.
142 143

  return energy_excess;
144

145
}