Simulate clusters of a track

[zhong]

Hello！
I want to get the clusters information of a track by using trackHeed, but it seems that Heed does not take the electric and magnetic field in the sensor into account for calculating the charged-particle trajectory.
I tried to use the EnableElectricField() &EnableMagneticField() function ,but it didn’t work.

If I want to simulate the relationship between the number of clusters and the electric and magnetic fields, what should I do

[jalopezg]

Hi @zhong; I am sure @hschindl can help you.

Cheers,
J.

[zhong]

OK，Thanks a lot！

[Crisps]

Can you show your code?

[hschindl]

In principle, Heed should take the electric/magnetic field into account when computing the trajectory but it might need a bit of fiddling with the stepping parameters.

Can you provide a minimal working example?

[zhong]

Here’s my working code.

#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <iomanip>
#include <cmath>
#include <vector>

#include <TCanvas.h>
#include <TMath.h>
#include <TROOT.h>
#include <TTree.h>
#include <TApplication.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TBuffer3D.h>
#include <TBuffer3DTypes.h>
#include <TVirtualViewer3D.h>
#include <TFile.h>
#include <TApplication.h>

#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/Plotting.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/TrackSrim.hh"

#include "Garfield/SolidBox.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include <time.h>

using namespace Garfield;
using namespace std;
using namespace ROOT;

int main() {
    bool debugging = true;
    double trackMom = 100.;
    int nTrk = 100;

    //save cluster information
    double xc = 0., yc = 0., zc = 0., tc = 0.;
    int nCluster;
    double trkion;
    //int Clus_posi;


    TFile *file = new TFile("ar_90_co2_10_0T_1GeV.root","RECREATE");
    TTree* c_position = new TTree("c","c Position");
    c_position->Branch("zc",&zc,"zc/D");
    c_position->Branch("ncl",&nCluster,"nCluster/I");
    c_position->Branch("trkion",&trkion,"trkion/D");
    
    //set gas information
    MediumMagboltz* gas = new MediumMagboltz();

    gas->SetComposition("ar",95.,"isobutane",2.,"CF4",3.);  
    gas->SetTemperature(293.15);
    gas->SetPressure(760.);
    gas->SetMaxElectronEnergy(200.);
    //gas->EnablePenningTransfer(0.44, 0.0, "He");
    gas->Initialise(true);

    // geo
    double zlength = 50.;
    SolidBox* box = new SolidBox(0., 0., 0., 20., 20, zlength); // center and half length
    GeometrySimple* geo = new GeometrySimple();
    geo->AddSolid(box, gas);

    ComponentAnalyticField* comp = new ComponentAnalyticField();
    if(debugging) comp->EnableDebugging();
    else comp->DisableDebugging();
    comp->SetGeometry(geo);
    comp->AddPlaneY(-20, 0., "pad_plane");
    comp->AddPlaneY(20, 10000., "HV");
    comp->SetMagneticField(0, 1, 0);
    

    Sensor* sensor = new Sensor();
    if(debugging) sensor->EnableDebugging();
    sensor->AddComponent(comp);
    sensor->AddElectrode(comp,"pad_plane");

    //setup a track 

    double x0 = 0, y0 =0 , z0 = -1*zlength, t0 = 0.; 
    //double dx0 = 0., dy0 = 1., dz0 = dy0/tan(theta);
    double dx0 = 0., dy0 = 0., dz0 = 1.;
    //double dy0 = sin(3.1415926*3/180) , dx0 = 0, dz0 = cos(3.1415926*3/180);
    TrackHeed* track = new TrackHeed();

    track->SetSensor(sensor);
    track->SetParticle("muon");
    track->SetMomentum(trackMom*1.e9);
    // track->EnableElectricField();
    // track->EnableMagneticField();
    // double maxrange = 0., rforstraight = 0., stepstraight = 0., stepcurved = 0.;
    // track->GetSteppingLimits(maxrange, rforstraight, stepstraight, stepcurved);
    // stepcurved = 0.02;
    // track->SetSteppingLimits(maxrange, rforstraight, stepstraight, stepcurved);


    for(int iTrk=0; iTrk<nTrk; iTrk++){
        track->NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);
        // Cluster coordinates
        //double xc = 0., yc = 0., zc = 0., tc = 0.;
        // Number of electrons produced in a collision
        int nc = 0;
        // Energy loss in a collision
        double ec = 0.;
        // Dummy variable (not used at present)
        double extra = 0.;
        //reset cluster
        nCluster = 0;
        //loop clusters 
        while(track->GetCluster(xc, yc, zc, tc, nc, ec, extra)){
            if((zc>= -50.) && (zc< 50.)){
            nCluster++;
            //Clus_posi = zc;
            c_position->Fill();
	        //cout << "cluster position : "<< zc << " deposited energy : "<< ec << endl;
            }
            
        }
        //cout<<"Cluster number: " << nCluster/100. <<endl;
        //trkion = nCluster/100.0;
        
    }
    c_position->Write(); 
    delete track;
    cout<<"end of the program"<<endl;
   
    file->Close();
    
}

[zhong]

yes, I have shown it, looking forward to your help

[hschindl]

At the momentum (100 GeV/c) and magnetic field (1 T) you are using, you wouldn’t expect a significant curvature of the track over the length of your detector (you can use the formula
p [GeV/c] = 0.3 B [T] r [m] to calculate the curvature radius).
If you lower the momentum you should be able to see a curved track (see the attached example).

test.C (1.3 KB)

[zhong]

Alright, Thanks for your help.