HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoAnalyseTrackPairs.cxx
Go to the documentation of this file.
1 //
2 // File HarpoAnalyseTrackPairs.cxx
3 //
11 #include "HarpoAnalyseTrackPairs.h"
12 #include "HarpoConfig.h"
13 #include "HarpoDetSet.h"
14 #include "HarpoDebug.h"
15 #include "HarpoDccEvent.h"
16 #include "Pmm2Event.h"
17 #include "HarpoEvent.h"
18 #include "HarpoRecoEvent.h"
19 #include "HarpoSimEvent.h"
20 
21 #include "TFile.h"
22 #include "TStyle.h"
23 #include "TCanvas.h"
24 #include "TLatex.h"
25 #include "TGraph.h"
26 #include "TF1.h"
27 #include "TMath.h"
28 #include "TLinearFitter.h"
29 #include "TSystem.h"
30 
31 #include <cstdlib>
32 #include <cstring>
33 #include <cassert>
34 #include <fstream>
35 #include <iostream>
36 
37 ClassImp(HarpoAnalyseTrackPairs);
38 
39 void HarpoAnalyseTrackPairs::print()
40  {
41  unsigned int nd; // number of detectors
42  HarpoEventHeader *hdr = fEvt->GetHeader();
43 
44  assert(hdr != NULL);
45  hdr->print();
46 
47  for (nd = 0; nd < gkNDetectors; nd++) {
48  // if (fEvt->isdataExist(nd)) {
49  HarpoDccMap *plane = fEvt->GetDccMap(nd);
50  if (plane != NULL )plane->print();
51  }
52  }
53 
54 void HarpoAnalyseTrackPairs::process()
55 {
56  // Bool_t drawEvent = kFALSE;
57  nEvents++;
58  if(nEvents%1000==0) Info("process","Event %li",nEvents);
59 
60 
61 
62  // if (gHDetSet->isExist(SIMDET)){
63  // HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
64  // Int_t nIonTr = simEvt->GetNtracks();
65  // TH1F* hTmp = new TH1F("hTmp","",512,0,512);
66  // for(Int_t i = 0; i<nIonTr; i++){
67  // TpcSimIonisationTrack* tr = simEvt->GetIonisationTrack(i);
68  // Int_t nPoints = tr->GetNpoints();
69  // for(Int_t j = 0; j<nPoints; j++){
70  // TpcSimIonisationPoint* point = tr->GetPoint(j);
71  // hTmp->Fill(point->GetZ(),point->GetNe());
72  // }
73  // }
74  // hTmp->Delete();
75  // }
76 
77 
78 
79 
80 
81 
82  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
83  // TClonesArray* clArray = reco->GetClustersArray();
84  // Int_t nCl = clArray->GetEntries();
85  TClonesArray* matchArray = reco->GetMatchingArray();
86  Int_t nMatch = matchArray->GetEntries();
87  TClonesArray* trArray = reco->GetTracksArray();
88  Int_t nTr = trArray->GetEntries();
89 
90  for(Int_t i = 0; i<nTr; i++){
91  HarpoRecoTracks* tr1 = (HarpoRecoTracks*)trArray->At(i);
92  Double_t xStart1 = tr1->GetXstart();
93  Double_t zStart1 = tr1->GetZstart();
94  Double_t xEnd1 = tr1->GetXend();
95  Double_t zEnd1 = tr1->GetZend();
96  Int_t plane1 = tr1->GetPlane();
97  Double_t dist = 1000000;
98  for(Int_t j = i+1; j<nTr; j++){
99  HarpoRecoTracks* tr2 = (HarpoRecoTracks*)trArray->At(j);
100  Int_t plane2 = tr2->GetPlane();
101  if(plane2 != plane1) continue;
102  Double_t xStart2 = tr2->GetXstart();
103  Double_t zStart2 = tr2->GetZstart();
104  Double_t xEnd2 = tr2->GetXend();
105  Double_t zEnd2 = tr2->GetZend();
106 
107  Double_t d[4];
108  d[0] = (xStart1-xStart2)*(xStart1-xStart2) + (zStart1-zStart2)*(zStart1-zStart2);
109  d[1] = (xEnd1-xEnd2)*(xEnd1-xEnd2) + (zEnd1-zEnd2)*(zEnd1-zEnd2);
110  d[2] = (xEnd1-xStart2)*(xEnd1-xStart2) + (zEnd1-zStart2)*(zEnd1-zStart2);
111  d[3] = (xStart1-xEnd2)*(xStart1-xEnd2) + (zStart1-zEnd2)*(zStart1-zEnd2);
112  for(Int_t k = 0; k<4; k++)
113  if(d[k]<dist) dist = d[k];
114 
115  hDistVertex->Fill(TMath::Sqrt(dist));
116  }
117  hDistVertex->Fill(400);
118  }
119 
120 
121  return;
122 
123  Int_t n3Dtr = 0;
124  for(Int_t i = 0; i<nMatch; i++){
125  HarpoRecoMatching* match = (HarpoRecoMatching*)matchArray->At(i);
126 
127  Int_t trX = match->GetiMatchTrX();
128  Int_t trY = match->GetiMatchTrY();
129  Int_t trX2 = match->GetiMatchTrX2();
130  Int_t trY2 = match->GetiMatchTrY2();
131 
132  n3Dtr++;
133  if(trX == trX2 && trY != trY2) n3Dtr++;
134  if(trX != trX2 && trY == trY2) n3Dtr++;
135  }
136 
137 
138  const Int_t kNdata = 18;
139  Double_t data[kNdata];
140  for(Int_t i = 0; i<kNdata; i++)
141  data[i] = -1000;
142  data[0] = reco->GetXstart();
143  data[1] = reco->GetYstart();
144  data[2] = reco->GetTstart();
145  data[12] = n3Dtr;
146 
147 
148 
149 
150  if (gHDetSet->isExist(SIMDET)){
151 
152  HarpoSimEvent *simEvt = (HarpoSimEvent *)(fEvt->GetDetEvent(SIMDET));
153  TpcSimMCEvent* mcevent = simEvt->GetMCEvent();
154  Int_t nMCtracks = mcevent->GetNtracks();
155  if(nMCtracks>=2){
156  for(Int_t i = 0; i<nMCtracks; i++){
157  TpcSimMCTrack* tr = mcevent->GetTrack(i);
158  Double_t p1[3];
159  p1[0] = tr->GetPx();
160  p1[1] = tr->GetPy();
161  p1[2] = tr->GetPz();
162  Double_t norm1 = 0;
163  for(Int_t k = 0; k<3; k++) norm1 += p1[k]*p1[k];
164  norm1 = TMath::Sqrt(norm1);
165  for(Int_t k = 0; k<3; k++) p1[k] /= norm1;
166  for(Int_t j = i+1; j<nMCtracks; j++){
167  TpcSimMCTrack* tr2 = mcevent->GetTrack(j);
168  Double_t p2[3];
169  p2[0] = tr2->GetPx();
170  p2[1] = tr2->GetPy();
171  p2[2] = tr2->GetPz();
172  Double_t norm2 = 0;
173  for(Int_t k = 0; k<3; k++) norm2 += p2[k]*p2[k];
174  norm2 = TMath::Sqrt(norm2);
175  for(Int_t k = 0; k<3; k++) p2[k] /= norm2;
176  Double_t u1u2 = 0;//, u1dX = 0, u2dX = 0;
177  for(Int_t k = 0; k<3; k++){
178  u1u2 += p1[k]*p2[k];
179  }
180  Double_t angle_open=acos(u1u2/(norm1*norm2));
181  Double_t angle_omega = atan2((p2[2]*p1[0]-p1[2]*p2[0]),(p2[2]*p1[1]-p1[2]*p2[1]));
182  hOmegaSim->Fill(angle_omega);
183  data[13] = angle_omega;
184  data[14] = angle_open;
185  }
186  }
187  }
188  }
189 
190 
191 
192 
193 
194  hN3Dtr->Fill(n3Dtr);
195  if(n3Dtr<2){
196  if(n3Dtr == 1){
197 
198  TClonesArray* trArrayK = reco->GetKalmanTracksArray();
199  Int_t nTrK = trArrayK->GetEntries();
200 
201  // Double_t x0tr[2][nTrK], chi2[2][nTrK];
202  Double_t slopetr[2][nTrK];
203  // Double_t xXtr[nTrK][2], xYtr[nTrK][2], zXtr[nTrK][2], zYtr[nTrK][2];
204  for(Int_t i = 0; i<nTrK; i++){
205  HarpoRecoKalmanTracks* track = (HarpoRecoKalmanTracks*)trArrayK->At(i);
206  Int_t idTr = track->GetNtrack();
207  Int_t plane = track->GetPlane();
208  if(plane!=XPLANE && plane !=YPLANE) continue;
209  if(idTr<0) continue;
210  if(idTr>=nTrK) continue;
211  // if(plane == XPLANE){
212  // xXtr[idTr][0] = track->GetXstart();
213  // zXtr[idTr][0] = track->GetZstart();
214  // xXtr[idTr][1] = track->GetXend();
215  // zXtr[idTr][1] = track->GetZend();
216  // }
217  // if(plane == YPLANE){
218  // xYtr[idTr][0] = track->GetXstart();
219  // zYtr[idTr][0] = track->GetZstart();
220  // xYtr[idTr][1] = track->GetXend();
221  // zYtr[idTr][1] = track->GetZend();
222  // }
223  // x0tr[plane][idTr] = track->GetX0();
224  slopetr[plane][idTr] = track->GetAngleTrack();
225  // chi2[plane][idTr] = track->GetChi2();
226  }
227 
228  for(Int_t i = 0; i<nMatch; i++){
229  HarpoRecoMatching* match = (HarpoRecoMatching*)matchArray->At(i);
230 
231  Int_t trX = match->GetiMatchTrX();
232  Int_t trY = match->GetiMatchTrY();
233  // Int_t trX2 = match->GetiMatchTrX2();
234  // Int_t trY2 = match->GetiMatchTrY2();
235 
236  // Double_t cosa = TMath::Cos(slopetr[0][trX])*TMath::Cos(slopetr[1][trY]);
237  // Double_t tanp = TMath::Tan(slopetr[0][trX])/TMath::Tan(slopetr[1][trY]);
238 
239 
240  // Double_t tanp = slopetr[0][trX]/slopetr[1][trY];
241  Double_t norm = slopetr[0][trX]*slopetr[0][trX]+slopetr[1][trY]*slopetr[1][trY] + 1;
242  Double_t cosa = TMath::Sqrt((norm-1.)/norm);
243 
244  data[9] = -10000;
245  data[10] = TMath::ATan2(slopetr[0][trX],slopetr[1][trY]);
246  data[11] = TMath::ACos(cosa);
247  }
248 
249  }
250  fNtuple->Fill(data);
251  return;
252  }
253  // Info("process","nMatch = %d",nMatch);
254  // if(nMatch!=2) return;
255  // TClonesArray* trArray = reco->GetTracksArray();
256  // Int_t nTr = trArray->GetEntries();
257  // if(nTr != 4) return;
258 
259 
260  TClonesArray* trArrayK = reco->GetKalmanTracksArray();
261  Int_t nTrK = trArrayK->GetEntries();
262 
263  Double_t x0tr[2][nTrK], slopetr[2][nTrK];//, chi2[2][nTrK];
264  Double_t xXtr[nTrK][2], xYtr[nTrK][2], zXtr[nTrK][2], zYtr[nTrK][2];
265  for(Int_t i = 0; i<nTrK; i++){
266  HarpoRecoKalmanTracks* track = (HarpoRecoKalmanTracks*)trArrayK->At(i);
267  Int_t idTr = track->GetNtrack();
268  Int_t plane = track->GetPlane();
269  if(plane!=XPLANE && plane !=YPLANE) continue;
270  if(idTr<0) continue;
271  if(idTr>=nTrK) continue;
272  for(Int_t j = i+1; j<nTrK; j++){
273  HarpoRecoKalmanTracks* track2 = (HarpoRecoKalmanTracks*)trArrayK->At(j);
274  Int_t plane2 = track2->GetPlane();
275  if(plane2!=plane) continue;
276  Int_t idTr2 = track2->GetNtrack();
277  if(idTr2<0) continue;
278  if(idTr2>=nTrK) continue;
279  Double_t dmin = 5000;
280  for(Int_t i1 = 0; i1<2; i1++){
281  Double_t x1, z1;
282  if(i1){
283  x1 = track->GetXstart();
284  z1 = track->GetZstart();
285  }else{
286  x1 = track->GetXend();
287  z1 = track->GetZend();
288  }
289  for(Int_t i2 = 0; i2<2; i2++){
290  Double_t x2, z2;
291  if(i2){
292  x2 = track->GetXstart();
293  z2 = track->GetZstart();
294  }else{
295  x2 = track->GetXend();
296  z2 = track->GetZend();
297  }
298  Double_t dX = TMath::Abs(x1 - x2);
299  Double_t dZ = TMath::Abs(z1 - z2);
300  if(dX + dZ < dmin){
301  dmin = dX + dZ;
302  // i1min = i1;
303  // i2min = i2;
304  }
305  }
306  }
307  hDends->Fill(dmin);
308  }
309  if(plane == XPLANE){
310  xXtr[idTr][0] = track->GetXstart();
311  zXtr[idTr][0] = track->GetZstart();
312  xXtr[idTr][1] = track->GetXend();
313  zXtr[idTr][1] = track->GetZend();
314  }
315  if(plane == YPLANE){
316  xYtr[idTr][0] = track->GetXstart();
317  zYtr[idTr][0] = track->GetZstart();
318  xYtr[idTr][1] = track->GetXend();
319  zYtr[idTr][1] = track->GetZend();
320  }
321  x0tr[plane][idTr] = track->GetX0();
322  slopetr[plane][idTr] = track->GetAngleTrack();
323  // chi2[plane][idTr] = track->GetChi2();
324  }
325 
326 
327  Int_t used[NTRACK][NTRACK];
328  for(Int_t i = 0; i<NTRACK; i++)
329  for(Int_t j = 0; j<NTRACK; j++)
330  used[i][j] = kBlack;
331  Int_t color = 2;
332  Double_t x0[nTrK], y0[nTrK], slopeX[nTrK], slopeY[nTrK];
333  Double_t xX[nTrK][2], xY[nTrK][2], zX[nTrK][2], zY[nTrK][2];
334  Double_t xEnds[nTrK][2], yEnds[nTrK][2], zEnds[nTrK][2];
335 
336  n3Dtr = 0;
337  for(Int_t i = 0; i<nMatch; i++){
338  HarpoRecoMatching* match = (HarpoRecoMatching*)matchArray->At(i);
339 
340  Int_t trX = match->GetiMatchTrX();
341  Int_t trY = match->GetiMatchTrY();
342  Int_t trX2 = match->GetiMatchTrX2();
343  Int_t trY2 = match->GetiMatchTrY2();
344  if(used[trX][trY] == kBlack){
345  used[trX][trY] = color;
346  color++;
347  }
348 
349 
350  for(Int_t k = 0; k<3; k++){
351  // Int_t iTrX, iTrY;
352  if(k == 1 && !(trX2 != trX && trY2 == trY)) continue;
353  if(k == 2 && !(trX2 == trX && trY2 != trY)) continue;
354  // switch(k){
355  // case 0: iTrX = trX; iTrY = trY; break;
356  // case 1: iTrX = trX2; iTrY = trY; break;
357  // case 2: iTrX = trX; iTrY = trY2; break;
358  // }
359  slopeX[n3Dtr] = slopetr[0][trX];
360  xX[n3Dtr][0] = xXtr[trX][0];
361  zX[n3Dtr][0] = zXtr[trX][0];
362  xX[n3Dtr][1] = xXtr[trX][1];
363  zX[n3Dtr][1] = zXtr[trX][1];
364  x0[n3Dtr] = x0tr[0][trX];
365 
366  slopeY[n3Dtr] = slopetr[1][trY];
367  xY[n3Dtr][0] = xYtr[trY][0];
368  zY[n3Dtr][0] = zYtr[trY][0];
369  xY[n3Dtr][1] = xYtr[trY][1];
370  zY[n3Dtr][1] = zYtr[trY][1];
371  y0[n3Dtr] = x0tr[1][trY];
372 
373  Double_t zDiffMin = 100000;
374  Int_t iMin = -1, jMin = -1;
375  iMin = 0;
376  if(zX[n3Dtr][1]<zX[n3Dtr][0]) iMin = 1;
377  jMin = 0;
378  if(zY[n3Dtr][1]<zY[n3Dtr][0]) jMin = 1;
379 
380  zDiffMin = TMath::Abs(zX[n3Dtr][iMin] - zY[n3Dtr][jMin]);
381  hDiffZstart->Fill(zDiffMin);
382  xEnds[n3Dtr][0] = xX[n3Dtr][iMin];
383  yEnds[n3Dtr][0] = xY[n3Dtr][jMin];
384  zEnds[n3Dtr][0] = 0.5*(zX[n3Dtr][iMin] + zY[n3Dtr][jMin]);
385  xEnds[n3Dtr][1] = xX[n3Dtr][1-iMin];
386  yEnds[n3Dtr][1] = xY[n3Dtr][1-jMin];
387  zEnds[n3Dtr][1] = 0.5*(zX[n3Dtr][1-iMin] + zY[n3Dtr][1-jMin]);
388 
389 
390  hStartFit->Fill(zEnds[n3Dtr][0],x0[n3Dtr] + slopeX[n3Dtr]*zEnds[n3Dtr][0],y0[n3Dtr] + slopeY[n3Dtr]*zEnds[n3Dtr][0]);
391 
392 
393  n3Dtr++;
394  }
395  }
396 
397  if(gHarpoDebug>1)
398  Info("process","N3dTracks = %d",n3Dtr);
399 
400  data[12] = n3Dtr;
401  Double_t distPocaVertexMin = 1e10, omega = 1000, open = 1000;
402 
403 
404  for(Int_t i = 0; i<n3Dtr; i++){
405  Double_t p1[3] = {slopeX[i], slopeY[i], 1};
406  Double_t x1[3] = {x0[i], y0[i], 0};
407  Double_t norm1 = 0;
408  for(Int_t k = 0; k<3; k++) norm1 += p1[k]*p1[k];
409  norm1 = TMath::Sqrt(norm1);
410  for(Int_t k = 0; k<3; k++) p1[k] /= norm1;
411  for(Int_t j = i+1; j<n3Dtr; j++){
412  Double_t p2[3] = {slopeX[j], slopeY[j], 1};
413  Double_t x2[3] = {x0[j], y0[j], 0};
414  Double_t norm2 = 0;
415  for(Int_t k = 0; k<3; k++) norm2 += p2[k]*p2[k];
416  norm2 = TMath::Sqrt(norm2);
417  for(Int_t k = 0; k<3; k++) p2[k] /= norm2;
418  Double_t u1u2 = 0, u1dX = 0, u2dX = 0;
419  for(Int_t k = 0; k<3; k++){
420  u1u2 += p1[k]*p2[k];
421  u1dX += p1[k]*(x1[k]-x2[k]);
422  u2dX += p2[k]*(x1[k]-x2[k]);
423  }
424 
425  Double_t t1 = (u1u2*u2dX - u1dX)/(1 - u1u2*u1u2);
426  Double_t t2 = - (u1u2*u1dX - u2dX)/(1 - u1u2*u1u2);
427 
428 
429  Double_t xPoca1[3], xPoca2[3], xPoca[3];
430  Double_t doca = 0;
431  for(Int_t k = 0; k<3; k++){
432  xPoca1[k] = x1[k] + t1*p1[k];
433  xPoca2[k] = x2[k] + t2*p2[k];
434  xPoca[k] = 0.5*(xPoca1[k] + xPoca2[k]);
435  doca += (xPoca1[k]-xPoca2[k])*(xPoca1[k]-xPoca2[k]);
436  }
437  doca = TMath::Sqrt(doca);
438  if(doca>0.01){
439  hDoca->Fill(doca);
440  }
441 
442 
443 
444  Double_t distMin = 100000;
445  Int_t iMin = -1, jMin = -1;
446  for(Int_t iTmp = 0; iTmp<2; iTmp++){
447  for(Int_t jTmp = 0; jTmp<2; jTmp++){
448  // Double_t distTmp = TMath::Abs(xEnds[i][iTmp] - xEnds[j][jTmp]) + TMath::Abs(yEnds[i][iTmp] - yEnds[j][jTmp]) + TMath::Abs(zEnds[i][iTmp] - zEnds[j][jTmp]);
449  Double_t distTmp = (xEnds[i][iTmp] - xEnds[j][jTmp])*(xEnds[i][iTmp] - xEnds[j][jTmp]) + (yEnds[i][iTmp] - yEnds[j][jTmp])*(yEnds[i][iTmp] - yEnds[j][jTmp]) + (zEnds[i][iTmp] - zEnds[j][jTmp])*(zEnds[i][iTmp] - zEnds[j][jTmp]);
450  if(distTmp<distMin){
451  distMin = distTmp;
452  iMin = iTmp;
453  jMin = jTmp;
454  }
455  }
456  }
457  hDistVertex->Fill(TMath::Sqrt(distMin));
458  // Double_t distPocaVertexI = TMath::Abs(xEnds[i][iMin] - xPoca[0]) + TMath::Abs(yEnds[i][iMin] - xPoca[1]) + TMath::Abs(zEnds[i][iMin] - xPoca[2]);
459  Double_t distPocaVertexI = (xEnds[i][iMin] - xPoca[0])*(xEnds[i][iMin] - xPoca[0]) + (yEnds[i][iMin] - xPoca[1])*(yEnds[i][iMin] - xPoca[1]) + (zEnds[i][iMin] - xPoca[2])*(zEnds[i][iMin] - xPoca[2]);
460  hDistPocaVertex->Fill(TMath::Sqrt(distPocaVertexI));
461  // Double_t distPocaVertexJ = TMath::Abs(xEnds[j][iMin] - xPoca[0]) + TMath::Abs(yEnds[j][iMin] - xPoca[1]) + TMath::Abs(zEnds[j][iMin] - xPoca[2]);
462  Double_t distPocaVertexJ = (xEnds[j][jMin] - xPoca[0])*(xEnds[j][jMin] - xPoca[0]) + (yEnds[j][jMin] - xPoca[1])*(yEnds[j][jMin] - xPoca[1]) + (zEnds[j][jMin] - xPoca[2])*(zEnds[j][jMin] - xPoca[2]);
463  hDistPocaVertex->Fill(TMath::Sqrt(distPocaVertexJ));
464 
465 
466 
467  hPoca->Fill(xPoca[2],xPoca[0],xPoca[1]);
468 
469 
470  Double_t angle_open=acos(u1u2/(norm1*norm2));
471  Double_t angle_omega = atan2((p2[2]*p1[0]-p1[2]*p2[0]),(p2[2]*p1[1]-p1[2]*p2[1]));
472  // hOpenAngle->Fill(angle_open*180./TMath::Pi());
473  // hOmegaAngle->Fill(angle_omega*180./TMath::Pi());
474 
475  Double_t distPocaVertex = (0.5*(xEnds[i][iMin] + xEnds[j][jMin]) - xPoca[0])*(0.5*(xEnds[i][iMin] + xEnds[j][jMin]) - xPoca[0]) + (0.5*(yEnds[i][iMin] + yEnds[j][jMin]) - xPoca[1])*(0.5*(yEnds[i][iMin] + yEnds[j][jMin]) - xPoca[1]) + (0.5*(zEnds[i][iMin] + zEnds[j][jMin]) - xPoca[2])*(0.5*(zEnds[i][iMin] + zEnds[j][jMin]) - xPoca[2]);
476  //hOmegaVsDistPocaVertex->Fill(angle_omega*180./TMath::Pi(),TMath::Sqrt(distPocaVertex));
477 
478  if(distPocaVertex<distPocaVertexMin){
479  distPocaVertexMin = distPocaVertex;
480  omega = angle_omega;
481  open = angle_open;
482  data[3] = xEnds[i][iMin];
483  data[4] = yEnds[i][iMin];
484  data[5] = zEnds[i][iMin];
485  data[6] = xEnds[j][jMin];
486  data[7] = yEnds[j][jMin];
487  data[8] = zEnds[j][jMin];
488  data[9] = distPocaVertex;
489  data[10] = omega;
490  data[11] = open;
491  for(Int_t k = 0; k<3; k++)
492  data[15+k] = 0.5*(p1[k]+p2[k]);
493  }
494 
495  }
496  }
497  hOmegaVsDistPocaVertex->Fill(omega*180./TMath::Pi(),TMath::Sqrt(distPocaVertexMin));
498  if(distPocaVertexMin<4900){
499  hOpenAngle->Fill(open*180./TMath::Pi());
500  hOmegaAngle->Fill(omega*180./TMath::Pi());
501  }
502 
503 
504 
505  hOmegaSimReco->Fill(data[13],data[10]);
506 
507 
508 
509 
510  fNtuple->Fill(data);
511 
512  // if(nEvents%1000==0)
513  // Save();
514 
515 }
516 
517 void HarpoAnalyseTrackPairs::Init()
518 {
519 
520  Double_t pi = TMath::Pi();
521  hStartFit = new TH3F("hStartFit",";z;x;y",200,0,511,100,0,288,100,0,288);
522  hPoca = new TH3F("hPoca",";z;x;y",200,0,511,100,0,288,100,0,288);
523  hDoca = new TH1F("hDoca",";DOCA [mm]",500,0,50);
524  hOpenAngle = new TH1F("hOpenAngle",";#theta [^{o}]",360,-180,180);
525  hOmegaAngle = new TH1F("hOmegaAngle",";#omega [^{o}]",360,-180,180);
526  hOmegaSim = new TH1F("hOmegaSim",";#omega [^{o}]",360,-pi,pi);
527  hOmegaSimReco = new TH2F("hOmegaSimReco",";#omega_{SIM} [^{o}];#omega_{RECO} [^{o}]",360,-pi,pi,360,-pi,pi);
528  hOmegaVsDistPocaVertex = new TH2F("hOmegaVsDistPocaVertex",";#theta [^{o}];dist(POCA,start track)",360,-180,180,500,0,500);
529  hDiffZstart = new TH1F("hDiffZstart",";z_{start,X} - z_{start,Y}",500,0,500);
530  hDistVertex = new TH1F("hDistVertex",";dist(start track X, start track Y)",500,0,500);
531  hDistPocaVertex = new TH1F("hDistPocaVertex",";dist(POCA,start track)",500,0,500);
532  hN3Dtr = new TH1F("hN3Dtr",";number of 3D tracks",50,0,50);
533  hDends = new TH1F("hDends",";number of 3D tracks",500,0,100);
534 
535  fNtuple = new TNtupleD("fNtuple","Vertex data","x0:y0:z0:xs1:ys1:zs1:xs2:ys2:zs2:dpv:omega:alpha:n3Dtr:omegasim:alphasim:px:py:pz");
536  fNtuple->SetDirectory(0);
537 
538  //
539  // fTree = new TNtuple("vertexTree","tree with vertex information","poca:doca:sXX:sZX:eXX:eZX:sXY:sZY:
540 
541 }
542 
543 void HarpoAnalyseTrackPairs::Save(char * /* mode */)
544  {
545 
546  OpenHistFile("trackpairs");
547 
548  hStartFit->Write();
549  hPoca->Write();
550  hDoca->Write();
551  hOpenAngle->Write();
552  hOmegaAngle->Write();
553  hDiffZstart->Write();
554  hDistVertex->Write();
555  hDistPocaVertex->Write();
556  hOmegaVsDistPocaVertex->Write();
557  hN3Dtr->Write();
558  hDends->Write();
559  hOmegaSim->Write();
560  hOmegaSimReco->Write();
561 
562  fNtuple->Write();
563 
564  }
565 
NTRACK
#define NTRACK
Definition: HarpoHoughTracking.h:22
HarpoEvent::GetRecoEvent
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
HarpoRecoEvent::GetKalmanTracksArray
TClonesArray * GetKalmanTracksArray()
Definition: HarpoRecoEvent.h:596
XPLANE
#define XPLANE
Definition: HarpoConfig.h:24
HarpoAnalyseTrackPairs::print
void print()
Overloaded method which do all job.
Definition: HarpoAnalyseTrackPairs.cxx:39
TpcSimMCTrack::GetPz
Double_t GetPz()
Definition: TpcSimMCEvent.h:32
HarpoRecoTracks::GetXend
Double_t GetXend()
Definition: HarpoRecoEvent.h:143
HarpoDetSet::isExist
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
Definition: HarpoDetSet.h:33
HarpoRecoTracks::GetAngleTrack
Double_t GetAngleTrack()
Definition: HarpoRecoEvent.h:136
HarpoRecoTracks
Track object, containing position, angle, charge and quality information.
Definition: HarpoRecoEvent.h:125
HarpoAnalyseTrackPairs
Dummy analysis to run as test and example. Give basic histograms of the data.
Definition: HarpoAnalyseTrackPairs.h:23
HarpoRecoEvent
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
Definition: HarpoRecoEvent.h:474
HarpoRecoMatching
Matching object, containing matched track number, and quality info.
Definition: HarpoRecoEvent.h:379
HarpoRecoEvent::GetMatchingArray
TClonesArray * GetMatchingArray()
Definition: HarpoRecoEvent.h:604
HarpoDccMap::print
virtual void print()
Definition: HarpoDccMap.cxx:124
TpcSimMCTrack::GetPy
Double_t GetPy()
Definition: TpcSimMCEvent.h:31
HarpoEventHeader
FullEvent Header not scecific to the detectors The class is ....
Definition: HarpoEventHeader.h:27
HarpoSimEvent
A class store HARPO row DCC event data and header. End provide access metods to the row data...
Definition: HarpoSimEvent.h:24
HarpoEventHeader::print
virtual void print()
Definition: HarpoEventHeader.cxx:21
HarpoRecoKalmanTracks
HarpoRecoTracks object, Obtained with Kalman filter.
Definition: HarpoRecoEvent.h:223
HarpoRecoTracks::GetX0
Double_t GetX0()
Definition: HarpoRecoEvent.h:137
HarpoAnalyse::OpenHistFile
TFile * OpenHistFile(const char *ananame)
Definition: HarpoAnalyse.cxx:150
HarpoDccMap
unpacked dcc data The class contains the data map for DCC or Feminos The data is stored as a 2D TMatr...
Definition: HarpoDccMap.h:29
HarpoEvent::GetDetEvent
HarpoDetEvent * GetDetEvent(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:93
SIMDET
Data from Keller temperuture and pressure sensors.
Definition: HarpoDet.h:22
HarpoRecoTracks::GetZstart
Double_t GetZstart()
Definition: HarpoRecoEvent.h:142
TpcSimMCEvent::GetNtracks
Int_t GetNtracks()
Definition: TpcSimMCEvent.h:100
TpcSimMCEvent::GetTrack
TpcSimMCTrack * GetTrack(Int_t i)
Definition: TpcSimMCEvent.h:91
gHarpoDebug
Long64_t gHarpoDebug
Definition: HarpoDebug.cxx:9
YPLANE
#define YPLANE
Definition: HarpoConfig.h:25
HarpoAnalyseTrackPairs::hDiffZstart
TH1F * hDiffZstart
Redefine empty default.
Definition: HarpoAnalyseTrackPairs.h:36
HarpoSimEvent::GetMCEvent
TpcSimMCEvent * GetMCEvent()
Definition: HarpoSimEvent.h:57
HarpoEvent::GetHeader
HarpoEventHeader * GetHeader()
Definition: HarpoEvent.cxx:80
HarpoAnalyse::fEvt
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
HarpoAnalyse::nEvents
ULong_t nEvents
Definition: HarpoAnalyse.h:75
HarpoRecoTracks::GetXstart
Double_t GetXstart()
Definition: HarpoRecoEvent.h:141
gkNDetectors
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
HarpoRecoEvent::GetTracksArray
TClonesArray * GetTracksArray()
Definition: HarpoRecoEvent.h:597
HarpoEvent::GetDccMap
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
TpcSimMCTrack::GetPx
Double_t GetPx()
Definition: TpcSimMCEvent.h:30
HarpoRecoTracks::GetZend
Double_t GetZend()
Definition: HarpoRecoEvent.h:144
gHDetSet
R__EXTERN HarpoDetSet * gHDetSet
Definition: HarpoDetSet.h:71
Generated by   doxygen 1.8.6