Particle direction

faca87 · May 31, 2020, 7:24am

Hello, yesterday, my supervisor wrote me to calculate the particle directions in an experimental setup. I replied asking me how I can do it and he replied

Take 2 points in x i.e. on two Silicon detecots on the same arm of the detector (for example the 35 and 37) putting the relative distance.

This is the experimental setup where you can see the axis directions, the subdetectors positions and the distances between them

Moreover, my supervisor’s colleague (he was in CC in the e.mail) also replied

To extrapolate, first you need the laboratory coordinate (because in the ROOT file there are just the loca ones).
For each event fill a vector using the hit coordinate in subdetectors 33, 35 and 37. For the error to assign to the hit (needed for the fit) use 0.030cm for all the trackers. If you get more hits whitin a tracker you or you have to reject by using for a reason, or you have to do all possible combinations and he sent me:

This macro for the coordinate traslation
TBAlign.C (13.5 KB)
TBAlign.h (4.3 KB)
This file for the fit.

fit_line.C (1.2 KB)

Unfortunately is the first time that I do a tracking reconstruction, so I’ve no idea about I have to do…I hope someone can help me to do it…

Anyway, I started to have a look to the macro, so

In TBAlign.h I repleaced

TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("si.root");
      if (!f || !f->IsOpen()) {
         f = new TFile("si.root");

by

 TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("C:/si-calo-sep18/si-500596.root");
      if (!f || !f->IsOpen()) {
         f = new TFile("C:/si-calo-sep18/si-500596.root");

In TBAlingn.c
I repleaced
TFile *newfile = new TFile("Tracker_info.root","recreate");
by
TFile *newfile = new TFile("C:/si-calo-sep18/si-500596-trasl.root","recreate");

but when I run the macro, I get

Does anyone know how why I get the error

(TBAlign) @0xc17e708

and how to do get the particle positions?
This is the ROOT file si-500596.root (1.4 MB)
Thank you

Please read tips for efficient and successful posting and posting code

ROOT Version: 6.20/24
Platform: Not Provided
Compiler: Not Provided

Wile_E_Coyote · May 31, 2020, 8:34am

That’s not any error (it’s an address which is returned).
You should have used: root [0] .L TBAlign.C++g

Actually, the “TBAlign.[hC]” files seem to come from the TTree::MakeClass and unfortunately some description has been cut out. Here’s the essential TBAlign::Loop() method:

void TBAlign::Loop()
{
//   In a ROOT session, you can do:
//      root> .L TBAlign.C
//      root> TBAlign t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
   }
}

BTW. All this should have been explained to you by the one who gave you these files.

faca87 · May 31, 2020, 9:27am

Thanks @Wile_E_Coyote

I runned the macro by using .L TBAlign.C++g and I don’t see anything on the terminal

moreover I don’t get the si-500596-trasl.root in the C:/si-calo-sep18/ directory.

About what you said TTree::MakeClass no…he didn’t wrote me about that and really I don’t know how to do the tracking reconstruction because it’s the first time… so please, if you know it, please can help me do to it? …
Thanks

Wile_E_Coyote · May 31, 2020, 9:31am

Execute commands shown in my previous post (given after “root>”).

faca87 · May 31, 2020, 9:43am

Ok. I executed the commands

root [0] .L TBAlign.c++g
root [1] TBAlign t
root [2] t.GetEntry(12);
root [3] t.Show();
root [4] t.Show(16);
root [5] t.Loop();

and this is all the log in ROOT

ROOT.txt (4.8 KB)

moreover I get this fit

and the new ROOT file

si-500596-trasl.root (1.7 MB)

Now, do you know what should I do?

Wile_E_Coyote · May 31, 2020, 9:54am

I guess, in the “TBAlign.C” file, in the “TBAlign::Loop” method, in its internal “for” loop over tree entries, you need to write some code which does what you want.

faca87 · May 31, 2020, 1:14pm

@Wile_E_Coyote, my supervisor just wrote me to calculate the particle direction in experimental setup and when I asked him how to do it, he just replied

Take 2 points in x i.e. on two Silicon detecots on the same arm of the detector (for example the 35 and 37) putting the relative distance .

So I think that I’ve to get the tracks of particles moving from the tracker 35 to the tracker 37 (shown in the image in first post).

Anyway, the macro just fill for the subdetectors 10,20,30,31 instead my supervisor and his collegue wrote me to study the trackers 33,35 and 37 so

In TBAlign.c I added the code
a.

TH1F *hnx33=new TH1F("hnx33","N hits in Si33, X coord",10,0,10);
  TH1F *hnx35=new TH1F("hnx35","N hits in Si35, X coord",10,0,10);
  TH1F *hnx37=new TH1F("hnx37","N hits in Si37, X coord",10,0,10);

b.
TH1F *hx33=new TH1F("hx33","X cood, Si33", 100,0,L33);
TH1F *hx35=new TH1F("hx35","X cood, Si35", 100,0,L35);
TH1F *hx37=new TH1F("hx37","X cood, Si37", 100,0,L37);

c.

TH2F *h2x33v31=new TH2F("h2x33v31","X, x33:x31", 100, 0, L31, 100, 0, L33);
  TH2F *h2x35v33=new TH2F("h2x35v33","X, x35:x33", 100, 0, L33, 100, 0, L35);
  TH2F *h2x37v35=new TH2F("h2x37v35","X, x37:x35", 100, 0, L35, 100, 0, L37);

d.

TH2F *h2y33v31=new TH2F("h2y33v31","Y, y33:y31", 100, 0, L31, 100, 0, L33);
   TH2F *h2y35v33=new TH2F("h2y35v33","Y, y35:y33", 100, 0, L33, 100, 0, L35);
   TH2F *h2y37v35=new TH2F("h2y37v35","Y, y37:y35", 100, 0, L35, 100, 0, L37);

e.

TH1F *h33diff = new TH1F("h33diff","X: x33-x31", 100, -1., 1.);
  TH1F *h35diff = new TH1F("h35diff","X: x35-x33", 100, -1., 1.);
  TH1F *h37diff = new TH1F("h37diff","X: x37-x35", 100, -1., 1.);

f.

TH1F *h33xz = new TH1F("h33xz","X: x33-x31, mrad", 100, -3., 3.);
    TH1F *h35xz = new TH1F("h35xz","X: x35-x33, mrad", 100, -3., 3.);
    TH1F *h37xz = new TH1F("h37xz","X: x37-x35, mrad", 100, -3., 3.);
TH1F *h33yz = new TH1F("h33yz","Y: y33-y31, mrad", 100, -3., 3.);
    TH1F *h35yz = new TH1F("h35yz","Y: y35-y33, mrad", 100, -3., 3.);
    TH1F *h37yz = new TH1F("h37yz","Y: y37-y35, mrad", 100, -3., 3.);

g.

TH2F *h2x33v31all=new TH2F("h2x33v31all","X, x33:x31", 100, 0, L31, 100, 0, L33);
     TH2F *h2x35v33all=new TH2F("h2x35v33all","X, x35:x33", 100, 0, L33, 100, 0, L35);
     TH2F *h2x37v35all=new TH2F("h2x37v35all","X, x37:x35", 100, 0, L35, 100, 0, L37);

h.

 TH2F *h2y33v31all=new TH2F("h2y33v31all","Y, y33:y31", 100, 0, L31, 100, 0, L33);
  TH2F *h2y35v33all=new TH2F("h2y35v33all","Y, y35:y33", 100, 0, L33, 100, 0, L35);
  TH2F *h2y37v35all=new TH2F("h2y37v35all","Y, y37:y35", 100, 0, L35, 100, 0, L37);

i.

hnx33->Fill(nx33);
       hnx35->Fill(nx35);
       hnx37->Fill(nx37);

l.

 if(nx31==1 && nx33==1) h2x33v31->Fill(x31[0],x33[0]);
   if(nx33==1 && nx35==1) h2x35v33->Fill(x33[0],x35[0]);
   if(nx35==1 && nx37==1) h2x37v35->Fill(x35[0],x37[0]);

m.

 h2x33v31all->Fill(x31[j],x33[j]);
             h2x35v33all->Fill(x33[j],x35[j]);
             h2x37v35all->Fill(x35[j],x37[j]);

m.

if(ny31==1 && ny33==1) h2y33v31->Fill(y31[0],y33[0]);
        if(ny33==1 && ny35==1) h2y35v33->Fill(y33[0],y35[0]);
        if(ny35==1 && ny37==1) h2y37v35->Fill(y35[0],y37[0]);

n.

h2y33v31all->Fill(y31[j],y33[j]);
           h2y35v33all->Fill(y33[j],y35[j]);
           h2y37v35all->Fill(y35[j],y37[j]);

o.

 if(nx31==1 && nx33==1) h33diff->Fill(x33[0]-x31[0]);
      if(nx33==1 && nx35==1) h35diff->Fill(x35[0]-x33[0]);
	  if(nx35==1 && nx37==1) h37diff->Fill(x37[0]-x35[0]);

p.
if(nx31==1 && nx33==1) h33xz->Fill(1e3*(x33[0]-x31[0])/(z33[0]-z31[0]));
if(nx33==1 && nx35==1) h35xz->Fill(1e3*(x35[0]-x33[0])/(z35[0]-z33[0]));
if(nx35==1 && nx37==1) h37xz->Fill(1e3*(x37[0]-x35[0])/(z37[0]-z35[0]));

   if(ny31==1 && ny33==1) h33yz->Fill(1e3*(y33[0]-y31[0])/(z33[0]-z31[0]));
    if(ny33==1 && ny35==1) h35yz->Fill(1e3*(y35[0]-y33[0])/(z35[0]-z33[0]));
     if(ny35==1 && ny37==1) h37yz->Fill(1e3*(y37[0]-y35[0])/(z37[0]-z35[0]));

q.

 h33xz->Fit("gaus");
   h33yz->Fit("gaus");
   h35xz->Fit("gaus");
   h35yz->Fit("gaus");
   h37xz->Fit("gaus");
   h37yz->Fit("gaus");

and this is the modified macro
TBAlign.C (19.8 KB)
TBAlign.h (4.3 KB)

I runned the macro getting this ROOT file si-500596-trasl.root (1.7 MB)

but it looks like there are problems…

In the 3 plots hx33, hx35, hx37 I read in the stat box that there are events, but I don’t see them, for example

image862×622 31.3 KB
In the 3 plots h2x33v31, h2x35v33, h2x37v35 I read in the stat box that there are events, but I don’t see them, for example

image865×611 33.3 KB
In the plot h33xz I read there are events, but I can’t see them

image856×622 28.9 KB
In the 3 plots h33diff, h35diff, h37diff, I read in the stat box that there are events, but I don’t see them, for example

Wile_E_Coyote · May 31, 2020, 1:25pm

Try to create your histograms with “automatic bins”, e.g.:

TH1F *hnx33 = new TH1F("hnx33", "N hits in Si33, X coord", 10, 0., 0.);

TH2F *h2x33v31 = new TH2F("h2x33v31", "X, x33:x31", 100, 0., 0., 100, 0., 0.);

faca87 · May 31, 2020, 2:22pm

Thank you @Wile_E_Coyote, So

I replaced

TH1F *hx33=new TH1F("hx33","X cood, Si33", 100,0,L33);
  TH1F *hx35=new TH1F("hx35","X cood, Si35", 100,0,L35);
  TH1F *hx37=new TH1F("hx37","X cood, Si37", 100,0,L37);

by

 TH1F *hx33=new TH1F("hx33","X cood, Si33", 100,0.,0.);
  TH1F *hx35=new TH1F("hx35","X cood, Si35", 100,0.,0.);
  TH1F *hx37=new TH1F("hx37","X cood, Si37", 100,0.,0.);

I replaced

TH2F *h2x33v31=new TH2F("h2x33v31","X, x33:x31", 100, 0, L31, 100, 0, L33);
  TH2F *h2x35v33=new TH2F("h2x35v33","X, x35:x33", 100, 0, L33, 100, 0, L35);
  TH2F *h2x37v35=new TH2F("h2x37v35","X, x37:x35", 100, 0, L35, 100, 0, L37);

by
TH2F *h2x33v31=new TH2F(“h2x33v31”,“X, x33:x31”, 100, 0., 0., 100, 0., 0.);
TH2F *h2x35v33=new TH2F(“h2x35v33”,“X, x35:x33”, 100, 0., 0., 100, 0., 0.);
TH2F *h2x37v35=new TH2F(“h2x37v35”,“X, x37:x35”, 100, 0., 0., 100, 0., 0.);

I replaced

 TH1F *h33xz = new TH1F("h33xz","X: x33-x31, mrad", 100, -3., 3.);
    TH1F *h35xz = new TH1F("h35xz","X: x35-x33, mrad", 100, -3., 3.);
    TH1F *h37xz = new TH1F("h37xz","X: x37-x35, mrad", 100, -3., 3.);

by

 TH1F *h33xz = new TH1F("h33xz","X: x33-x31, mrad", 100, 0., 0.);
    TH1F *h35xz = new TH1F("h35xz","X: x35-x33, mrad", 100, 0., 0.);
    TH1F *h37xz = new TH1F("h37xz","X: x37-x35, mrad", 100, 0., 0.);

I replaced

TH1F *h33diff = new TH1F("h33diff","X: x33-x31", 100, -1., 1.);
  TH1F *h35diff = new TH1F("h35diff","X: x35-x33", 100, -1., 1.);
  TH1F *h37diff = new TH1F("h37diff","X: x37-x35", 100, -1., 1.);

by

TH1F *h33diff = new TH1F("h33diff","X: x33-x31", 100, 0., 0.);
  TH1F *h35diff = new TH1F("h35diff","X: x35-x33", 100, 0., 0.);
  TH1F *h37diff = new TH1F("h37diff","X: x37-x35", 100, 0., 0.);

and now I can see the graph.

Now I have the plots that my supervisor’s collegue talked about, i.e. the hit coordinates i silicon detectors 33, 35 and 37.

But isn’t strange that I get negative coordinates?

Anyway, Given that he also sent me this code fit_line.C (1.2 KB)
and he wrote me to assign error 0.030cm to the hit position for the fit, I guess that somehow I should add this code in macro ROOT and fit the coordinate…Do you know what I should do?

Wile_E_Coyote · May 31, 2020, 2:27pm

Show your current progress to your supervisor and his colleague and ask for further guidance / help.

faca87 · May 31, 2020, 2:36pm

Yes, @Wile_E_Coyote, you are right…it’s better first show them the plots.

But I was thinking that I used the code that you showed me

root> .L TBAlign.C
      root> TBAlign t
      root> t.GetEntry(12); // Fill t data members with entry number 12
      root> t.Show();       // Show values of entry 12
      root> t.Show(16);     // Read and show values of entry 16
      root> t.Loop();       // Loop on all entries

Writing t.GetEntry(12) and t.Show(16); did I miss something? I mean that I that I specified 12 and 16.

Thanks

Wile_E_Coyote · May 31, 2020, 2:38pm

“ROOT User’s Guide” -> “Trees” -> “Using TTree::MakeClass”

faca87 · May 31, 2020, 3:15pm

I tried to run changing the numbers I see that the x coordinate plots don’t change.

faca87 · June 1, 2020, 7:49am

@Wile_E_Coyote I wrote to my supervisor’s and his colleague showing what we did yesterday. My supervisor’s colleague replied that it’s right to get negative coordinates and he wrote me what I must do now. I’m writing you what he wrote in the email

You can do a fit by using the 3 trackers 33,35 and 37. This is what you have to do:

float x[3]; float z[3]; float ex[3];
ex[0]=0.1; ex[1]=0.1; ex[2]=0.1;

// where ex[] are the errors on each hit. Actually the errors are 0.03mm, but take 0.1mm so that you take into account also a not perfect alignment

Fill the x and z arrays by xh and zh coordinates

float a,b,ae,be,chi2,q;
      fit_line(z,x,3,ex,a,b,ae,be,chi2,q);

//The fit function is x = a + b*z
//ae and be are the errors releated to a,b parameters, we don't neet q.

That’s what he wrote me yesterday night. So, I added:

On the top, before the TBAlign::Loop()

#include "fit_line.C"
float x[3]; float z[3]; float ex[3];
   float a,b,ae,be,chi2,q;
   ex[0]=0.1; ex[1]=0.1; ex[2]=0.1;

Inside the for(Long64_t jentry=0; jentry<nentries;jentry++) {}
a. To fill the x and z arrays:

  	if(nx33==1){
     		x[0]=x33[0];
     		z[0]=z33[0];	
		 } 
     	if(nx35==1){
     		x[1]=x35[0];
     		z[1]=z35[0];
		 }
     	if(nx37==1) {
     		x[2]=x37[0];
     		z[2]=z37[0];
		 }

On the bottom, inside TBAlign::Loop() , but out of for(Long64_t jentry=0; jentry<nentries;jentry++) {}

fit_line(z,x,3,ex,a,b,ae,be,chi2,q);

but when I run I get this error

Wile_E_Coyote · June 1, 2020, 10:33am

Move the line which generates this error to inside of the TBAlign::Loop or try:
float x[3]; float z[3]; float ex[3] = {0.1, 0.1, 0.1};

faca87 · June 1, 2020, 12:21pm

Thank you @Wile_E_Coyote I moved the line

ex[0]=0.1; ex[1]=0.1; ex[2]=0.1;

into the TBAlign::Loop and I runned it looks like working (I get the ROOT file) , but during the run I read

Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x
Error in <TInterpreter::TCling::AutoLoad>: failure loading library TBAlign_c for x

is it important?

This is the full log log.txt (9.3 KB)

Moreover, how can see now the results (graph and parameters values) of the fit

fit_line(z,x,3,ex,a,b,ae,be,chi2,q);

?

Thank you

fit_line.C (1.2 KB)
TBAlign.C (21.6 KB)
TBAlign.h (4.3 KB)

Wile_E_Coyote · June 1, 2020, 12:55pm

@Axel I have no idea where these errors may come from.

What concerns “results”, you can, for example, create several std::vector and / or TGraphErrors and / or TGraph2DErrors and / or TTree and add entries to them after each call to fit_line (using the values returned in a,b,ae,be,chi2,q).

faca87 · June 1, 2020, 2:07pm

Thank you @Wile_E_Coyote, I tried the TGraphErrors

I included #include <TGraphErrors.h>
I added

fit_line(z,x,3,ex,a,b,ae,be,chi2,q);
    auto c0 = new TCanvas("c0","Graph Si33",200,10,700,500);
    auto gr0 = new TGraphErrors(0,x,z,ex);
  	gr0->SetTitle("Si33");
   	gr0->SetMarkerColor(4);
   	gr0->SetMarkerStyle(21);
   	gr0->Fit("fit_line");
   	gr0->Draw("E");
   	c0->Print("C:/si-calo-sep18/Si33.pdf");
   	
    auto c1 = new TCanvas("c1","Graph Si35",200,10,700,500);
    auto gr1 = new TGraphErrors(1,x,z,ex);
  	gr1->SetTitle("Si35");
   	gr1->SetMarkerColor(4);
   	gr1->SetMarkerStyle(21);
   	gr1->Fit("fit_line");
   	gr1->Draw("E");
   	c1->Print("C:/si-calo-sep18/Si35.pdf");
   	
   	 
    auto c2 = new TCanvas("c2","Graph Si37",200,10,700,500);
    auto gr2 = new TGraphErrors(2,x,z,ex);
  	gr2->SetTitle("Si37");
   	gr2->SetMarkerColor(4);
   	gr2->SetMarkerStyle(21);
   	gr2->Fit("fit_line");
   	gr2->Draw("E");
   	c2->Print("C:/si-calo-sep18/Si37.pdf");

But I get 3 empty plots ( I mean a with canvas), for example
Si33.pdf (12.2 KB)

Wile_E_Coyote · June 1, 2020, 7:49pm

Try to play with this example macro:

{
  float z[] = {1873.93, 1995.43, 2125.43}; // z coordinates of 33, 35, 37
  float x[] = {-7.0, -13.5, -20.5}; // x coordinates of hits (an example)
  float ex[] = {0.1, 0.1, 0.1}; // errors of x coordinates of hits
  int n = 3; // number of data points (2 or 3)
  
  gStyle->SetOptFit(1112);
  TCanvas *c = new TCanvas("c", "track", 200, 10, 700, 500);
  TGraphErrors *g = new TGraphErrors(n, z, x, 0, ex); // hits x versus z
  g->SetTitle("track;z;x");
  g->SetMarkerColor(4);
  g->SetMarkerStyle(4);
  const char *f_name = "pol1"; // standard ROOT linear function
  g->Fit(f_name, ""); // e.g. "" or "Q" or "V"
  TF1 *f = g->GetFunction(f_name); // retrieve the fit function
  if (f) {
    std::cout << "fit function name = " << f->GetName() << std::endl;
    std::cout << "Chi2 = " << f->GetChisquare() << std::endl;
    std::cout << "NDF = " << f->GetNDF() << std::endl;
    std::cout << "Prob = " << f->GetProb() << std::endl;
    for (Int_t i = 0; i < f->GetNpar(); i++) {
      std::cout << f->GetParName(i) << " = " << f->GetParameter(i)
                << " +- " << f->GetParError(i) << std::endl;
    }
  }
  g->Draw("AP");
  c->SaveAs("track.pdf");
  // delete c; // no longer needed
  // delete g; // no longer needed
}

faca87 · June 2, 2020, 9:10am

Thank you @Wile_E_Coyote it worked!
I first fitted, by a 2 gaussian function, the x coordinates for the three silicon trackers, then, using your macro I fitted the mean values got by the 2 gaussian fit in function of the z values!
Just a question…

Given that the x errors is so small (0.1 cm), using the

g->Draw("AP");

I just see a square

so I tried the

1.g->Draw("ap");

but I get the same plot

g->Draw("[]"); optiion.
but I get an empty plot…

It isn’t a big problem…just to know if there is a way to get it better.

TBAlign.C (30.0 KB) TBAlign.h (4.3 KB)