C++ Sensor::GetElectronSignal方法代码示例

//.........这里部分代码省略.........


		int nElastic;
		int nIonising;
		int nAttachment;
		int nInelastic;
		int nExcitation;
		int nSuperelastic;
		gas-> GetNumberOfElectronCollisions(nElastic,nIonising,nAttachment,nInelastic,nExcitation,nSuperelastic);
		std::cout<<"nElastic="<<nElastic<<std::endl;
		std::cout<<"nIonising="<<nIonising<<std::endl;
		std::cout<<"nAttachment="<<nAttachment<<std::endl;
		std::cout<<"nInelastic="<<nInelastic<<std::endl;
		std::cout<<"nExcitation="<<nExcitation<<std::endl;
		std::cout<<"nSuperelastic="<<nSuperelastic<<std::endl;


		tot_nElastic += nElastic;
		tot_nIonising += nIonising;
		tot_nAttachment += nAttachment;
		tot_nInelastic += nInelastic;
		tot_nExcitation += nExcitation;
		tot_nSuperelastic += nSuperelastic;

		// To get Signal from each strip we get the it from each sensor
		// by integrating over time

		for (int j=0; j<n_strips_x; j++)
		{	
			const char* name =   ("Strip"+to_string(j)).c_str();

			for(int i = 0; i < nTimeBins; i++)
			{
				double value = sensor->GetElectronSignal(name, i);
				
				if (std::isnan(value) == false )
				{
					QstripsTot[j] += value;
					QtimeTot[i] += value;
				}
			}    
		}

		
		for(int i = 0; i < nTimeBins; i++)
		{
			QtimeIonTot[i] += sensor->GetIonSignal("StripIon", i);
		}
		
		//// Sanity check
		//double sum_q_time = 0;
		//double sum_q_strip = 0;
		
		//for (int j=0; j<n_strips_x; j++)
			//sum_q_strip += QstripsTot[j];
			
		//for(int i = 0; i < nTimeBins; i++)
			//sum_q_time  += QtimeTot[i];
			
		//if (sum_q_strip != sum_q_time)
		//{
			//cout << "WARNING: sum_q_strip != sum_q_time" << endl;
			//cout << "sum_q_strip: " << sum_q_strip << endl;
			//cout << "sum_q_time: "  << sum_q_time  << endl;
		//}

//.........这里部分代码省略.........
  aval->EnableMagneticField();
  // Additional optional switches
  //aval->EnableExcitationMarkers();
  //aval->EnableIonisationMarkers();


  //Add ionizing particle using Heed
  // Here we add a negative pion with some momentum, units in [eV/c]
  const double momentum = 300.e+06; // eV/c
  TrackHeed* track = new TrackHeed();
  track->SetParticle("pi-");
  track->SetMomentum(momentum);
  track->SetSensor(sensor);
  track->EnableMagneticField();
  track->EnableElectricField();
  track->EnablePlotting(viewdrift);
  
  // Cluster info
  double xcls, ycls, zcls, tcls, e, extra;
  xcls = ycls = zcls = tcls = e = extra = -999.;
  // Electron info
  double xele, yele, zele, tele, eele, dxele, dyele, dzele;
  // Electron start and endpoints, momentum direction and status
  double x0ele, y0ele, z0ele, t0ele, e0ele;// start point
  double x1ele, y1ele, z1ele, t1ele, e1ele;// endpoint
  double dx1ele, dy1ele, dz1ele; // momentum direction
  int status1ele; // status
  int n = 0; // number of electrons in cluster
  bool cls_ret;// return OK if cluster is OK

  // The initial impact position of the incoming ionising track
  double track_x = 0.15;// [cm]
  double track_y = 0.3;
  double track_z = 0.0;
  // Momentum direction of incoming track
  double track_dx = 0.0;
  double track_dy = -1.0; // Track coming downstream
  double track_dz = 0.0;

  // Now create a single track
  track->NewTrack(track_x, track_y, track_z, tMin, track_dx, track_dy, track_dz);
  cls_ret = track->GetCluster(xcls, ycls, zcls, tcls, n, e, extra);
  
  // To accumulate the electron signal for each strip separately
  double Qstrip1 = 0.0; double Qstrip2 = 0.0; double Qstrip3 = 0.0;

  // Now loop over electrons in cluster
  for(int j = 1; j <= n; j++){
    track->GetElectron(j-1, xele, yele, zele, tele, eele, dxele, dyele, dzele);

    // Simulate an avalanche for the current electron
    aval->AvalancheElectron(xele, yele, zele, tele, eele, dxele, dyele, dzele); 
  }

  // To get Signal from each strip we get the it from each sensor
  // by integrating over time
  for(int i = 0; i < nTimeBins; i++){
    Qstrip1 += fabs(sensor->GetElectronSignal("Strip1", i));
    Qstrip2 += fabs(sensor->GetElectronSignal("Strip2", i));
    Qstrip3 += fabs(sensor->GetElectronSignal("Strip3", i));
  }    

  // Now calculate the reconstructed pion position in the Micromegas
  // using the weighted average of the Strip signals
  double mean = 0.0;
  mean = (Qstrip1*Xstrip1 + Qstrip2*Xstrip2 + Qstrip3*Xstrip3)/(Qstrip1+Qstrip2+Qstrip3);
  std::cout << "---------------------------" << std::endl;
  std::cout << "XMean: " << mean << " cm, XTrue: " << track_x << " cm" << std::endl;
  std::cout << "XTrue - XMean: " << fabs(track_x - mean)*10000.0 << " um" << std::endl;
  std::cout << "---------------------------" << std::endl;

  // -- Plotting --

  // Now plot the drift lines
  viewdrift->Plot();

  // View the Field
  ViewField * viewfield = new ViewField();
  viewfield->SetComponent(cmpAmp);
  viewfield->SetSensor(sensor);
  viewfield->SetCanvas((TCanvas*)c->cd(2));
  viewfield->SetWeightingFieldRange(0.0, 10000.0);
  c->cd(2);
  viewfield->PlotContour();


  // View the strip signals
  TCanvas * c2 = new TCanvas("c2", "c2", 1000, 20, 700, 500);
  ViewSignal* signalView = new ViewSignal();
  signalView->SetSensor(sensor);
  signalView->SetCanvas(c2);
  if(Qstrip1 > 0)signalView->PlotSignal("Strip1");
  if(Qstrip2 > 0)signalView->PlotSignal("Strip2");
  if(Qstrip3 > 0)signalView->PlotSignal("Strip3");


  app.Run();
  return 0;

}