FEDRA emulsion software from the OPERA Collaboration
tr_carbon_new.C File Reference

Functions

void BookHistsV ()
 
void ClearHistsV ()
 
void d ()
 
void DrawHistsV ()
 
void DrawHlist (TCanvas *c, TObjArray h)
 
void FillHistsV (Vertex &v)
 
void td ()
 
void tr_carbon_new (const char *rcut="nseg>2&&t.eFlag>-1")
 
void vd (int trmin=3, int trmax=100, float amin=.015, float amax=2.)
 

Variables

float dpp = .5
 
EdbDataProcdproc =0
 
EdbDisplayds =0
 
EdbPVRecgAli =0
 
EdbVertexRecgEVR =0
 
TObjArray hld1
 
TObjArray hld2
 
TObjArray hld3
 
TObjArray hld4
 
TObjArray hlist
 
TH1F * hp [23] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}
 
int MINMULT = 3
 
float momentum = 1.
 
int nsegMin = 3
 
float ProbMinP = 0.01
 
float ProbMinV = 0.01
 
bool ReFit = false
 
bool ReProp = false
 
bool usemom = false
 
bool UseSegPar = false
 

Function Documentation

◆ BookHistsV()

void BookHistsV ( )

198 {
199  if (!hp[3]) hp[3] = new TH1F("Hist_Vt_Sx","Vertex X sigma",100,0.,10.);
200  if (!hp[4]) hp[4] = new TH1F("Hist_Vt_Sy","Vertex Y sigma",100,0.,10.);
201  if (!hp[5]) hp[5] = new TH1F("Hist_Vt_Sz","Vertex Z sigma",100,0.,500.);
202  if (!hp[6]) hp[6] = new TH1F("Hist_Vt_Chi2","Vertex Chi-square/D.F.",25,0.,5.);
203  if (!hp[7]) hp[7] = new TH1F("Hist_Vt_Prob","Vertex Chi-square Probability",26,0.,1.04);
204  if (!hp[9]) hp[9] = new TH1F("Hist_Vt_Nvert","Number of reconstructed vertexs",100,0.,100.);
205  if (!hp[22]) hp[22] = new TH1F("Hist_Vt_Ntracks","Number of tracks at vertex", 20, 0.,20.);
206  if (!hp[14]) hp[14] = new TH1F("Hist_Vt_Dist","RMS track-vertex distance",100, 0.,50.);
207  if (!hp[18]) hp[18] = new TH1F("Hist_Vt_AngleV","RMS track-track angle, mrad", 100, 0.,1000.);
208 
209  if (!hp[10]) hp[10] = new TH1F("Hist_Vt_Ntrack","Number of reconstructed tracks",100,0.,10000.);
210  if (!hp[11]) hp[11] = new TH1F("Hist_Vt_Nsegs","Number of track segments",60,0.,60.);
211  if (!hp[12]) hp[12] = new TH1F("Hist_Vt_Dtrack","Track DeltaX, microns",100,-5.,+5.);
212  if (!hp[17]) hp[17] = new TH1F("Hist_Vt_Trprob","Track probability",110, 0.,1.1);
213  if (!hp[20]) hp[20] = new TH1F("Hist_Vt_Angle","Track angle, mrad", 100, 0.,1000.);
214  if (!hp[21]) hp[21] = new TH1F("Hist_Vt_Momen","Track momentum, GeV/c", 100, 0.,20.);
215 
216  hld1.Add(hp[3]);
217  hld1.Add(hp[4]);
218  hld1.Add(hp[5]);
219  hld1.Add(hp[6]);
220  hld1.Add(hp[7]);
221  hld1.Add(hp[9]);
222  hld1.Add(hp[22]);
223  hld1.Add(hp[14]);
224  hld1.Add(hp[18]);
225 
226  hld2.Add(hp[10]);
227  hld2.Add(hp[11]);
228  hld2.Add(hp[12]);
229  hld2.Add(hp[17]);
230  hld2.Add(hp[20]);
231  hld2.Add(hp[21]);
232 
233  for (int i=0; i<22; i++)
234  if (hp[i]) hlist.Add(hp[i]);
235 }
TH1F * hp[23]
Definition: tr_carbon_new.C:22
TObjArray hld1
Definition: tr_carbon_new.C:23
TObjArray hld2
Definition: tr_carbon_new.C:23
TObjArray hlist
Definition: tr_carbon_new.C:23

◆ ClearHistsV()

void ClearHistsV ( )

337 {
338  for(int i=0; i<22; i++) {
339  if (hp[i]) hp[i]->Clear();
340  }
341 }

◆ d()

void d ( )

334 { DrawHistsV();}
void DrawHistsV()
Definition: tr_carbon_new.C:250

◆ DrawHistsV()

void DrawHistsV ( )

251 {
252 
253  int n = hld1.GetEntries();
254  if (n)
255  {
256  TCanvas *cv1 = new TCanvas("Vertex_reconstruction_1","MC Vertex reconstruction", 760, 760);
257  DrawHlist(cv1, hld1);
258  }
259 
260  n = hld2.GetEntries();
261  if (n)
262  {
263  TCanvas *cv2 = new TCanvas("Vertex_reconstruction_2","Vertex reconstruction (tracks hists)", 600, 760);
264  DrawHlist(cv2, hld2);
265  }
266 
267  n = hld3.GetEntries();
268  if (n)
269  {
270  TCanvas *cv3 = new TCanvas("Vertex_reconstruction_3","Vertex reconstruction (eff hists)", 600, 760);
271  DrawHlist(cv3, hld3);
272  }
273 
274  n = hld4.GetEntries();
275  if (n)
276  {
277  TCanvas *cv4 = new TCanvas("Vertex_reconstruction_4","Vertex reconstruction (ntracks)", 600, 760);
278  DrawHlist(cv4, hld4);
279  }
280 }
TCanvas * cv4
Definition: RecDispMC_Profiles.C:63
TCanvas * cv1
Definition: RecDispMC_Profiles.C:63
TCanvas * cv2
Definition: RecDispMC_Profiles.C:63
TCanvas * cv3
Definition: RecDispMC_Profiles_new.C:70
new TCanvas()
void DrawHlist(TCanvas *c, TObjArray h)
Definition: tr_carbon_new.C:282
TObjArray hld4
Definition: tr_carbon_new.C:23
TObjArray hld3
Definition: tr_carbon_new.C:23

◆ DrawHlist()

void DrawHlist ( TCanvas c,
TObjArray  h 
)

283 {
284  int n = h.GetEntries();
285  if (n < 2)
286  {
287  }
288  else if (n < 3)
289  {
290  c->Divide(1,2);
291  }
292  else if (n < 4)
293  {
294  c->Divide(1,3);
295  }
296  else if (n < 5)
297  {
298  c->Divide(2,2);
299  }
300  else if (n < 6)
301  {
302  c->Divide(2,3);
303  }
304  else if (n < 7)
305  {
306  c->Divide(2,3);
307  }
308  else if (n < 8)
309  {
310  c->Divide(2,4);
311  }
312  else if (n < 9)
313  {
314  c->Divide(2,4);
315  }
316  else if (n < 10)
317  {
318  c->Divide(3,3);
319  }
320  else if (n < 11)
321  {
322  c->Divide(2,5);
323  }
324  else
325  {
326  c->Divide(3,5);
327  }
328  for(int i=0; i<n; i++) {
329  c->cd(i+1);
330  if (h.At(i)) h.At(i)->Draw();
331  }
332 }

◆ FillHistsV()

void FillHistsV ( Vertex v)

239 {
240  hp[3]->Fill(v.vxerr());
241  hp[4]->Fill(v.vyerr());
242  hp[5]->Fill(v.vzerr());
243  hp[6]->Fill(v.ndf() > 0 ? v.chi2()/v.ndf() : 0);
244  hp[7]->Fill(v.prob());
245  hp[14]->Fill(v.rmsDistAngle());
246  hp[18]->Fill(v.angle()*1000.);
247  hp[22]->Fill(v.ntracks());
248 }
double vyerr() const
$\sqrt{\sigma_{vy}^2}$ vertex $y$-error
unsigned short int ndf() const
degrees of freedom of vertex fit
Definition: VtVertex.C:238
float chi2() const
$\chi^2$ of vertex fit
Definition: VtVertex.C:236
double vzerr() const
$\sqrt{\sigma_{vz}^2}$ vertex $z$-error
double angle() const
unsigned short int ntracks() const
no of tracks in vertex $n$
Definition: VtVertex.C:240
double rmsDistAngle() const
calc rms dist and rms angle
Definition: VtVertex.C:1073
float prob() const
upper tail $\chi^2$ probability
Definition: VtVertex.C:237
double vxerr() const
$\sqrt{\sigma_{vx}^2}$ vertex $x$-error

◆ td()

void td ( )
146 {
147  TObjArray *trarr=dset->PVR()->eTracks;
148  gStyle->SetPalette(1);
149 
150  const char *title = "display-tracks";
151  if(!(ds=EdbDisplay::EdbDisplayExist(title)))
152  ds=new EdbDisplay(title,-50000.,50000.,-50000., 50000., -4000.,80000.);
153 
154  ds->SetVerRec(gEVR);
155  ds->SetDrawTracks(4);
156  ds->SetArrTr( trarr );
157  ds->Draw();
158 }
EdbDataProc * dset
Definition: RecDispEX.C:9
EdbPVRec * PVR() const
Definition: EdbDataSet.h:197
virtual void Draw(Option_t *option="")
Definition: EdbDisplayBase.cxx:786
Definition: EdbDisplay.h:21
void SetArrTr(TObjArray *arr)
Definition: EdbDisplay.cxx:430
void SetDrawTracks(int opt)
Definition: EdbDisplay.h:97
void SetVerRec(EdbVertexRec *evr)
Definition: EdbDisplay.h:91
static EdbDisplay * EdbDisplayExist(const char *title)
Definition: EdbDisplay.cxx:232
TObjArray * eTracks
Definition: EdbPVRec.h:161
EdbDisplay * ds
Definition: tr_carbon_new.C:5
EdbVertexRec * gEVR
Definition: tr_carbon_new.C:6

◆ tr_carbon_new()

void tr_carbon_new ( const char *  rcut = "nseg>2&&t.eFlag>-1")
27 {
28  BookHistsV();
29 
30  dproc = new EdbDataProc("Carbonium_data_set.def");
31  dproc->InitVolume(100, rcut);
32 
33  gAli = dproc->PVR();
34 
35  int ntr = gAli->eTracks->GetEntries();
36  printf("\n%d tracks are in the list\n",ntr);
37 
38  EdbTrackP *tra = 0;
39  int ntrgood = 0;
40  for(int i=0; i<ntr; i++)
41  {
42  tra = ((EdbTrackP*)(gAli->eTracks->At(i)));
43  if (tra->Flag() >= 0) ntrgood++;
44  }
45  printf("\n%d good tracks are in the list\n",ntrgood);
46 
47  if (ReFit)
48  {
49  for(int i=0; i<ntr; i++)
50  {
51  tra = ((EdbTrackP*)(gAli->eTracks->At(i)));
52  if (tra->Flag() < 0) continue;
56  }
57  gAli->FitTracks(momentum,0.139);
58  }
59 
60  gAli->FillCell(50,50,0.015,0.015);
61 
62  if (ReProp)
63  {
65  ntrgood = 0;
66  for(int i=0; i<ntr; i++)
67  {
68  tra = ((EdbTrackP*)(gAli->eTracks->At(i)));
69  if (tra->Flag() >= 0) ntrgood++;
70  }
71  printf("\n%d good tracks are in the list after additional propagation\n",
72  ntrgood);
73  }
74 
75  if (hp[10]) hp[10]->Fill(ntr);
76  float ang = 0.;
77  float tx = 0.;
78  float ty = 0.;
79  // tracks loop
80  ntr = gAli->eTracks->GetEntries();
81  int nt = 0;
82  for(int i=0; i<ntr; i++)
83  {
84  tra = (EdbTrackP*)(gAli->eTracks->At(i));
85  if (!tra) continue;
86 // tra->Print();
87  if (tra->Flag() != -10 && tra->N() >= 2) // skip "broken" tracks
88  {
89  // fill hist with number of segments
90  if (hp[11]) hp[11]->Fill(tra->N());
91  // fill hist with difference between measured and reconstructed coordinate
92  if (hp[12]) hp[12]->Fill((tra->GetSegment(0))->X()-(tra->GetSegmentF(0))->X());
93  // fill hist with track probability
94  if (hp[17]) hp[17]->Fill(tra->Prob());
95  tx = tra->TX();
96  ty = tra->TY();
97  ang = TMath::ACos(1./(1.+tx*tx+ty*ty))*1000.;
98  // fill hist with track Teta angle
99  if (hp[20]) hp[20]->Fill(ang);
100  // fill hist with track momentum
101  if (hp[21]) hp[21]->Fill(tra->P());
102  nt++;
103  }
104  }
105  printf("%d good tracks was found\n",nt);
106 
107  gEVR = new EdbVertexRec();
108  if (!(gAli->eTracks)) gAli->eTracks = new TObjArray(100);
110  if (!(gAli->eVTX)) gAli->eVTX = new TObjArray(100);
111  gEVR->eVTX = gAli->eVTX;
112  gEVR->ePVR = gAli;
113 
114  gEVR->eZbin = 100.; // default is 100. microns
115  gEVR->eAbin = 0.01; // default is 0.01 rad
116  gEVR->eDZmax = 3000.; // default is 3000. microns
117  gEVR->eProbMin = ProbMinV; // default 0.01, i.e 1%
118  gEVR->eImpMax = 55.; // default is 25. microns
119  gEVR->eUseSegPar = UseSegPar; // default is FALSE - use fitted track parameters
120  gEVR->eQualityMode= 0; // 0 - prob/(sx**2+sy**2), 1 = 1./aver_impact
121 
122  int nvtx = gEVR->FindVertex();
123  printf("%d 2-track vertexes was found\n",nvtx);
124 
125  if(nvtx != 0) {
126  int nadd = gEVR->ProbVertexN();
127  printf("%d vertexes with number of tracks >2 was found\n",nadd);
128  printf("----------------------------------------------\n");
129  EdbVertex *v = 0;
130  Vertex *vt = 0;
131  for(int i=0; i<nvtx; i++)
132  {
133  v = (EdbVertex*)(gEVR->eVTX->At(i));
134  if (!v) continue;
135 // if (i<10) v->Print();
136  if (v->Flag() < 0) continue;
137  vt = v->V();
138  if (vt) FillHistsV(*vt);
139  }
140  }
141 // vd();
142 }
Definition: EdbDataSet.h:180
int InitVolume(int datatype=0, const char *rcut="1")
Definition: EdbDataSet.cxx:2071
TObjArray * eVTX
Definition: EdbPVRec.h:162
void FillCell(float stepx, float stepy, float steptx, float stepty)
Definition: EdbPVRec.cxx:1088
int PropagateTracks(int nplmax, int nplmin, float probMin=0.05, int ngapMax=3, int design=0)
Definition: EdbPVRec.cxx:2483
void FitTracks(float p=10., float mass=0.139, TObjArray *gener=0, int design=0)
Definition: EdbPVRec.cxx:1889
Float_t Prob() const
Definition: EdbSegP.h:153
Float_t TX() const
Definition: EdbSegP.h:172
Float_t P() const
Definition: EdbSegP.h:149
Float_t TY() const
Definition: EdbSegP.h:173
void SetErrorP(float sp2)
Definition: EdbSegP.h:93
Int_t Flag() const
Definition: EdbSegP.h:146
Definition: EdbPattern.h:118
EdbSegP * GetSegmentFirst() const
Definition: EdbPattern.h:194
EdbSegP * GetSegmentLast() const
Definition: EdbPattern.h:195
Int_t N() const
Definition: EdbPattern.h:182
EdbSegP * GetSegmentF(int i) const
Definition: EdbPattern.h:201
EdbSegP * GetSegment(int i) const
Definition: EdbPattern.h:200
Int_t eQualityMode
Definition: EdbVertex.h:183
Float_t eAbin
Definition: EdbVertex.h:176
Bool_t eUseSegPar
Definition: EdbVertex.h:182
Float_t eImpMax
Definition: EdbVertex.h:179
Float_t eZbin
Definition: EdbVertex.h:175
Float_t eDZmax
Definition: EdbVertex.h:177
Float_t eProbMin
Definition: EdbVertex.h:178
Definition: EdbVertex.h:194
Int_t ProbVertexN()
Definition: EdbVertex.cxx:1448
TObjArray * eVTX
Definition: EdbVertex.h:205
Int_t FindVertex()
Definition: EdbVertex.cxx:1087
EdbPVRec * ePVR
Definition: EdbVertex.h:206
TObjArray * eEdbTracks
Definition: EdbVertex.h:204
Definition: EdbVertex.h:68
VERTEX::Vertex * V() const
Definition: EdbVertex.h:154
Int_t Flag() const
Definition: EdbVertex.h:124
Definition: VtVertex.hh:88
Double_t X
Definition: tlg2pattern.C:77
bool ReFit
Definition: tr_carbon_new.C:16
float momentum
Definition: tr_carbon_new.C:13
float ProbMinV
Definition: tr_carbon_new.C:8
void FillHistsV(Vertex &v)
Definition: tr_carbon_new.C:238
bool UseSegPar
Definition: tr_carbon_new.C:15
EdbDataProc * dproc
Definition: tr_carbon_new.C:3
float dpp
Definition: tr_carbon_new.C:14
bool ReProp
Definition: tr_carbon_new.C:18
void BookHistsV()
Definition: tr_carbon_new.C:197
float ProbMinP
Definition: tr_carbon_new.C:9
EdbPVRec * gAli
Definition: tr_carbon_new.C:4

◆ vd()

void vd ( int  trmin = 3,
int  trmax = 100,
float  amin = .015,
float  amax = 2. 
)
162 {
163  TObjArray *varr = new TObjArray();
164  TObjArray *tarr = new TObjArray();
165 
166  EdbVertex *v=0;
167  EdbTrackP *t=0;
168 
169  int nv = gEVR->eVTX->GetEntries();
170  for(int i=0; i<nv; i++) {
171  v = (EdbVertex *)(gEVR->eVTX->At(i));
172 
173  if( v->N()<trmin ) continue;
174  if( v->N()>trmax ) continue;
175  if( v->Flag()<0 ) continue;
176  if( v->N()<3 && v->MaxAperture()<amin ) continue;//
177  if( v->N()<3 && v->MaxAperture()>amax ) continue;
178 
179  varr->Add(v);
180  for(int j=0; j<v->N(); j++) tarr->Add( v->GetTrack(j) );
181  }
182 
183  gStyle->SetPalette(1);
184  const char *title = "display-vertices";
185  if(!(ds=EdbDisplay::EdbDisplayExist(title)))
186 // ds=new EdbDisplay(title,-30000.,2000., -2000.,42000.,-80000.,80000.);
187  ds=new EdbDisplay(title,60000.,70000., 45000.,50000., 0., 110000.);
188 
189  ds->SetVerRec(gEVR);
190  ds->SetArrTr( tarr );
191  ds->SetArrV( varr );
192  ds->SetDrawTracks(4);
193  ds->SetDrawVertex(1);
194  ds->Draw();
195 }
TTree * t
Definition: check_shower.C:4
void SetDrawVertex(int opt)
Definition: EdbDisplay.h:108
void SetArrV(TObjArray *arrv)
Definition: EdbDisplay.cxx:500
EdbTrackP * GetTrack(int i)
Definition: EdbVertex.h:141
Int_t N() const
Definition: EdbVertex.h:121
Float_t MaxAperture()
Definition: EdbVertex.cxx:252

Variable Documentation

◆ dpp

float dpp = .5

◆ dproc

EdbDataProc* dproc =0

◆ ds

EdbDisplay* ds =0

◆ gAli

EdbPVRec* gAli =0

◆ gEVR

EdbVertexRec* gEVR =0

◆ hld1

TObjArray hld1

◆ hld2

TObjArray hld2

◆ hld3

TObjArray hld3

◆ hld4

TObjArray hld4

◆ hlist

TObjArray hlist

◆ hp

TH1F* hp[23] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}

◆ MINMULT

int MINMULT = 3

◆ momentum

float momentum = 1.

◆ nsegMin

int nsegMin = 3

◆ ProbMinP

float ProbMinP = 0.01

◆ ProbMinV

float ProbMinV = 0.01

◆ ReFit

bool ReFit = false

◆ ReProp

bool ReProp = false

◆ usemom

bool usemom = false

◆ UseSegPar

bool UseSegPar = false