HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoMatchingVertex.cxx
Go to the documentation of this file.
1 //
2 // File HarpoMatchingVertex.cxx
3 //
9 #include "HarpoMatchingVertex.h"
10 #include "HarpoConfig.h"
11 #include "HarpoDetSet.h"
12 #include "HarpoDebug.h"
13 #include "HarpoDccEvent.h"
14 #include "Pmm2Event.h"
15 #include "HarpoEvent.h"
16 #include "HarpoRecoEvent.h"
17 #include "HarpoSimEvent.h"
18 #include "TpcSimMCEvent.h"
19 #include "MakeNiceHisto.h"
20 #include "HarpoTools.h"
21 
22 #include "TFile.h"
23 #include "TStyle.h"
24 #include "TCanvas.h"
25 #include "TRootEmbeddedCanvas.h"
26 #include "TGListBox.h"
27 #include "TLatex.h"
28 #include "TGraph.h"
29 #include "TF1.h"
30 #include "TMath.h"
31 #include "TGLabel.h"
32 #include "TRandom.h"
33 #include "TArrayD.h"
34 #include "TSystem.h"
35 
36 #include <cstdlib>
37 #include <cstring>
38 #include <cassert>
39 #include <fstream>
40 #include <iostream>
41 
42 ClassImp(HarpoMatchingVertex);
43 
45  {
46  unsigned int nd; // number of detectors
48 
49  assert(hdr != NULL);
50  hdr->print();
51 
52  for (nd = 0; nd < gkNDetectors; nd++) {
53  // if (fEvt->isdataExist(nd)) {
54  HarpoDccMap *plane = fEvt->GetDccMap(nd);
55  if (plane != NULL )plane->print();
56  }
57  }
58 
60 {
61  // Bool_t drawEvent = kFALSE;
62  nEvents++;
63 
64  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
65  if(!reco) return;
66  TClonesArray* vArray = reco->GetVertexArray();
67  if(!vArray) return;
68  Int_t nV = vArray->GetEntries();
69  if(nV<=0) return;
70  reco->ResetVertex3DArray();
71 
72 
73  HarpoDccMap* m[2];
74  m[0] = fEvt->GetDccMap(0);
75  m[1] = fEvt->GetDccMap(1);
76  // Initialise histograms here
77  for(Int_t i = 0; i<kMaxNvertex; i++){
78  for(Int_t j = 0; j<3; j++){ // profile for tracks 1 and 2. 0: both tracks
79  hQvT[i][j]->Reset();
80  }
81  hOverlap[i]->Reset();
82  }
83 
84  for(Int_t iV = 0; iV<nV; iV++){
85  HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV);
86  Int_t plane = vertex->GetPlane();
87  FillQvT(vertex, hQvT[iV][1], hQvT[iV][2], hOverlap[iV], m[plane]);
88  hQvT[iV][0]->Add(hQvT[iV][1],hQvT[iV][2]);
89  }
90 
91 
92  Double_t mean[7][nV*nV];
93  Double_t sigma[7][nV*nV];
94  Double_t cov[7][nV*nV];
95  Int_t nPair = 0;
96  Int_t iX[nV*nV], iY[nV*nV];
97  for(Int_t i = 0; i<nV*nV; i++){
98  iX[i] = -1;
99  iY[i] = -1;
100  for(Int_t j = 0; j<7; j++){
101  mean[j][i] = 0;
102  sigma[j][i] = 0;
103  cov[j][i] = -2;
104  }
105  }
106  for(Int_t iVx = 0; iVx<nV; iVx++){
107  HarpoRecoVertex* vertexX = (HarpoRecoVertex*)vArray->At(iVx);
108  if(vertexX->GetPlane()!=XPLANE) continue;
109  Double_t zVx = vertexX->GetZvertex();
110  for(Int_t iVy = 0; iVy<nV; iVy++){
111  HarpoRecoVertex* vertexY = (HarpoRecoVertex*)vArray->At(iVy);
112  if(vertexY->GetPlane()!=YPLANE) continue;
113  Double_t zVy = vertexY->GetZvertex();
114  if(TMath::Abs(zVx - zVy)>20) continue;
115 
116  cov[0][nPair] = CompareHist(hQvT[iVx][0],hQvT[iVy][0],hOverlap[iVx],hOverlap[iVy],mean[0][nPair],sigma[0][nPair]);
117 
118  if(vertexX->GetThetaVertex() != 0 && vertexY->GetThetaVertex() != 0){
119  cov[1][nPair] = CompareHist(hQvT[iVx][1],hQvT[iVy][1],hOverlap[iVx],hOverlap[iVy],mean[1][nPair],sigma[1][nPair]);
120  cov[2][nPair] = CompareHist(hQvT[iVx][2],hQvT[iVy][2],hOverlap[iVx],hOverlap[iVy],mean[2][nPair],sigma[2][nPair]);
121  cov[3][nPair] = CompareHist(hQvT[iVx][1],hQvT[iVy][2],hOverlap[iVx],hOverlap[iVy],mean[3][nPair],sigma[3][nPair]);
122  cov[4][nPair] = CompareHist(hQvT[iVx][2],hQvT[iVy][1],hOverlap[iVx],hOverlap[iVy],mean[4][nPair],sigma[4][nPair]);
123  cov[5][nPair] = Covariance(hQvT[iVx][1],hQvT[iVy][1],hQvT[iVx][2],hQvT[iVy][2],hOverlap[iVx],hOverlap[iVy]);
124  cov[6][nPair] = Covariance(hQvT[iVx][1],hQvT[iVy][2],hQvT[iVx][2],hQvT[iVy][1],hOverlap[iVx],hOverlap[iVy]);
125  // hCov->Fill(cov[1][nPair]);
126  // hCov->Fill(cov[2][nPair]);
127  // hCov->Fill(cov[3][nPair]);
128  // hCov->Fill(cov[4][nPair]);
129  //cov[5][nPair] = TMath::Max(cov[1][nPair],cov[2][nPair]);
130  //cov[6][nPair] = TMath::Max(cov[3][nPair],cov[4][nPair]);
131  hCov->Fill((cov[1][nPair]+cov[2][nPair])/(cov[3][nPair]+cov[4][nPair]));
132  hCov2->Fill(cov[1][nPair]+cov[2][nPair],cov[3][nPair]+cov[4][nPair]);
133  hCov3->Fill(cov[5][nPair],cov[6][nPair]);
134  }
135  iX[nPair] = iVx;
136  iY[nPair] = iVy;
137  nPair++;
138  }
139 
140  }
141 
142 
143  if(gHarpoDebug>1){
144  for(Int_t i = 0; i<nPair; i++){
145  Info("process","X: %d, Y:%d",iX[i],iY[i]);
146  for(Int_t j = 0; j<5; j++){
147  Info("process"," %d: %g %g",j,mean[j][i],sigma[j][i]);
148  }
149  }
150  }
151 
152  Int_t* index = new Int_t[nPair];
153  // TMath::Sort(nPair,sigma[0],index,kFALSE);
154  TMath::Sort(nPair,cov[0],index,kTRUE);
155 
156 
157 
158 
159 
160 
161  Double_t psim[2][3];
162  // Double_t angle_omega_sim = -1000;
163  if (gHDetSet->isExist(SIMDET)){
165  TpcSimMCEvent* mcEvent = simEvt->GetMCEvent();
166  Int_t nMCtr = mcEvent->GetNtracks();
167  Double_t norm[2] = {0,0};
168  for(Int_t i = 0; i<nMCtr; i++){
169  TpcSimMCTrack* track = mcEvent->GetTrack(i);
170  psim[i][0] = track->GetPy();
171  psim[i][1] = track->GetPx();
172  psim[i][2] = track->GetPz();
173  norm[i] = TMath::Sqrt(psim[i][0]*psim[i][0]+psim[i][1]*psim[i][1]+psim[i][2]*psim[i][2]);
174  //hPsim->Fill(track->GetPx(),track->GetPy(),track->GetPz());
175  }
176  // Double_t u1u2 = 0;
177  for(Int_t k = 0; k<3; k++){
178  psim[0][k] /= norm[0];
179  psim[1][k] /= norm[1];
180  }
181  //angle_omega_sim = atan2((psim[1][2]*psim[0][0]-psim[0][2]*psim[1][0]),(psim[1][2]*psim[0][1]-psim[0][2]*psim[1][1]));
182  //angle_omega_sim = 0.5*(TMath::ATan(
183 
184  }
185 
186 
187 
188 
189 
190 
191 
192 
193 
194 
195 
196 
197 
198 
199 
200 
201 
202 
203 
204 
205  Int_t used[nV];
206  for(Int_t i = 0; i< nV; i++) used[i] = 0;
207  for(Int_t i = 0; i<nPair; i++){
208  if(used[iX[index[i]]]) continue;
209  if(used[iY[index[i]]]) continue;
210  used[iX[index[i]]] = 1;
211  used[iY[index[i]]] = 1;
212 
213  if(gHarpoDebug>1)
214  Info("process","Keeping pair %d - %d",iX[index[i]],iY[index[i]]);
215 
216  Int_t j;
217  Double_t chi2 = 0;
218  Bool_t test = cov[5][index[i]] < cov[6][index[i]];
219  // Bool_t test = (cov[1][index[i]]+cov[2][index[i]] < cov[3][index[i]]+cov[4][index[i]]);
220  //Bool_t test = TMath::Max(cov[1][index[i]],cov[2][index[i]])<TMath::Max(cov[3][index[i]],cov[4][index[i]]);
221  if(test){
222  j = 1;
223  chi2 = cov[3][index[i]]+cov[4][index[i]];
224  if(gHarpoDebug>1)
225  Info("process","1->2 2->1");
226  }else{
227  if(gHarpoDebug>1)
228  Info("process","1->1 2->2");
229  j = 0;
230  chi2 = cov[1][index[i]]+cov[2][index[i]];
231  }
232  // if(cov[1][index[i]]+cov[2][index[i]] == cov[3][index[i]]+cov[4][index[i]])
233  if(cov[5][index[i]] == cov[6][index[i]])
234  j = Int_t(gRandom->Rndm()+0.5);
235 
236 
237 
238 
239 
240  HarpoRecoVertex* vertexX = (HarpoRecoVertex*)vArray->At(iX[index[i]]);
241  Double_t xVx = vertexX->GetXvertex();
242  Double_t zVx = vertexX->GetZvertex();
243  Double_t pxVx = vertexX->GetPxVertex();
244  Double_t pzVx = vertexX->GetPzVertex();
245  Double_t thVx = vertexX->GetThetaVertex();
246  Double_t cThX = TMath::Cos(thVx*0.5);
247  Double_t sThX = TMath::Sin(thVx*0.5);
248  Double_t pxX[2], pzX[2];
249  pxX[0] = pxVx*cThX + pzVx*sThX;
250  pzX[0] = pzVx*cThX - pxVx*sThX;
251  pxX[1] = pxVx*cThX - pzVx*sThX;
252  pzX[1] = pzVx*cThX + pxVx*sThX;
253  HarpoRecoVertex* vertexY = (HarpoRecoVertex*)vArray->At(iY[index[i]]);
254  Double_t xVy = vertexY->GetXvertex();
255  Double_t zVy = vertexY->GetZvertex();
256  Double_t pxVy = vertexY->GetPxVertex();
257  Double_t pzVy = vertexY->GetPzVertex();
258  Double_t thVy = vertexY->GetThetaVertex();
259  Double_t cThY = TMath::Cos(thVy*0.5);
260  Double_t sThY = TMath::Sin(thVy*0.5);
261  Double_t pxY[2], pzY[2];
262  pxY[0] = pxVy*cThY + pzVy*sThY;
263  pzY[0] = pzVy*cThY - pxVy*sThY;
264  pxY[1] = pxVy*cThY - pzVy*sThY;
265  pzY[1] = pzVy*cThY + pxVy*sThY;
266 
267  Double_t data[23];
268  data[21] = j;
269 
270  if(fPerfectMatching == 1 && gHDetSet->isExist(SIMDET)){
271  if( (psim[0][0]-psim[1][0])*(psim[0][1]-psim[1][1])*(pxX[0]-pxX[1])*(pxY[0]-pxY[1])>0 )
272  j = 0;
273  else
274  j = 1;
275  }
276  if(fPerfectMatching == 2)
277  j = Int_t(gRandom->Rndm()+0.5);
278 
279  Double_t p1[3];
280  if(pzX[0]>0 && pzY[j]>0){
281  p1[0] = pxX[0]/pzX[0];
282  p1[1] = pxY[j]/pzY[j];
283  p1[2] = 1;
284  }else if(pzX[0]<0 && pzY[j]<0){
285  p1[0] = -pxX[0]/pzX[0];
286  p1[1] = -pxY[j]/pzY[j];
287  p1[2] = -1;
288  }else{
289  p1[0] = 0;
290  p1[1] = 0;
291  p1[2] = 0;
292  }
293  Double_t p2[3];
294  if(pzX[1]>0 && pzY[1-j]>0){
295  p2[0] = pxX[1]/pzX[1];
296  p2[1] = pxY[1-j]/pzY[1-j];
297  p2[2] = 1;
298  }else if(pzX[1]<0 && pzY[1-j]<0){
299  p2[0] = -pxX[1]/pzX[1];
300  p2[1] = -pxY[1-j]/pzY[1-j];
301  p2[2] = -1;
302  }else{
303  p2[0] = 0;
304  p2[1] = 0;
305  p2[2] = 0;
306  }
307  Double_t norm1 = 0, norm2 = 0;
308  for(Int_t k = 0; k<3; k++){
309  norm1 += p1[k]*p1[k];
310  norm2 += p2[k]*p2[k];
311  }
312  norm1 = TMath::Sqrt(norm1);
313  norm2 = TMath::Sqrt(norm2);
314  for(Int_t k = 0; k<3; k++){
315  p1[k] /= norm1;
316  p2[k] /= norm2;
317  }
318 
319  data[0] = 0.5*(zVx+zVy);
320  data[1] = xVx;
321  data[2] = xVy;
322  data[3] = pxX[0];
323  data[4] = pxX[1];
324  data[5] = pxY[0];
325  data[6] = pxY[1];
326  data[7] = cov[0][index[i]];
327  data[8] = cov[1][index[i]];
328  data[9] = cov[2][index[i]];
329  data[10] = cov[3][index[i]];
330  data[11] = cov[4][index[i]];
331  data[12] = cov[5][index[i]];
332  data[13] = cov[6][index[i]];
333  data[14] = psim[0][0];
334  data[15] = psim[0][1];
335  data[16] = psim[0][2];
336  data[17] = psim[1][0];
337  data[18] = psim[1][1];
338  data[19] = psim[1][2];
339 
340  if( (psim[0][0]-psim[1][0])*(psim[0][1]-psim[1][1])*(pxX[0]-pxX[1])*(pxY[0]-pxY[1])>0 )
341  data[20] = 0;
342  else
343  data[20] = 1;
344  if((psim[0][0]-psim[1][0])*(psim[0][1]-psim[1][1])*(pxX[0]-pxX[1])*(pxY[0]-pxY[1])==0)
345  data[20] = -1;
346 
347  Double_t u1u2 = p1[0]*p2[0] + p1[1]*p2[1] + p1[2]*p2[2];
348  Double_t angle_open=acos(u1u2);
349  data[22] = angle_open;
350  fNtuple->Fill(data);
351 
352  HarpoRecoVertex3D* v3D = new HarpoRecoVertex3D(xVx, xVy, 0.5*(zVx+zVy), p1[0], p1[1], p1[2], p2[0], p2[1], p2[2], iX[index[i]], iY[index[i]], chi2);
353 
354  reco->AddVertex3D(v3D);
355  v3D->Delete();
356 
357  }
358 }
359 
360 Double_t HarpoMatchingVertex::CompareHist(TH1F* hX, TH1F* hY, TH1F* hoverlapX, TH1F* hoverlapY, Double_t &mean, Double_t &sigma)
361 {
362 
363  mean = -1;
364  sigma = -1;
365 
366  Int_t nbins = hX->GetXaxis()->GetNbins();
367  if(hY->GetXaxis()->GetNbins() != nbins) return -1;
368 
369 
370  TArrayD* arr = new TArrayD(nbins);
371  mean = 0;
372  sigma = 0;
373  Int_t n = 0;
374  Double_t W = 0, meanX = 0, sigX = 0, meanY = 0, sigY = 0, cov = 0;
375  Double_t qXold = -1000, qYold = -1000;
376  for(Int_t bin = 1; bin<=nbins; bin++){
377  Double_t qX = hX->GetBinContent(bin);
378  Double_t qY = hY->GetBinContent(bin);
379  Double_t dqX = -10000;
380  if(qXold>0) dqX = (qX - qXold);
381  Double_t dqY = -10000;
382  if(qYold>0) dqY = (qY - qYold);
383  if(fPerfectMatching==3){
384  dqX = qX;
385  dqY = qY;
386  }
387  if(fPerfectMatching==4){
388  if(qXold>0) dqX = (qX - qXold)/(qX+qXold);
389  if(qYold>0) dqY = (qY - qYold)/(qY+qYold);
390  }
391  qXold = qX;
392  qYold = qY;
393  if(qY<=0) continue;
394  if(qX<=0) continue;
395  if(dqY==-10000) continue;
396  if(dqX==-10000) continue;
397  Double_t w = 1;
398  if(hoverlapX->GetBinContent(bin)>0) w *= 0.1;
399  if(hoverlapY->GetBinContent(bin)>0) w *= 0.1;
400 
401 
402  arr->AddAt(qX/qY,bin);
403  //std::cout << qX/qY << " ";
404  mean += w*qX/qY;
405  sigma += w*qX/qY*qX/qY;
406 
407  // cov += qX*qY;
408  // meanX += qX;
409  // sigX += qX*qX;
410  // meanY += qY;
411  // sigY += qY*qY;
412 
413  cov += w*dqX*dqY;
414  meanX += w*dqX;
415  sigX += w*dqX*dqX;
416  meanY += w*dqY;
417  sigY += w*dqY*dqY;
418  W += w;
419  n++;
420  }
421 
422 
423 
424  // std::cout << std::endl;
425  if(W==0){
426  delete arr;
427  return 0;
428  }
429  meanX /= W;
430  sigX /= W;
431  sigX -= meanX*meanX;
432  meanY /= W;
433  sigY /= W;
434  sigY -= meanY*meanY;
435  cov /= W;
436  cov -= meanX*meanY;
437  cov /= TMath::Sqrt(sigX*sigY);
438 
439  mean /= W;
440  sigma = sigma/W - mean*mean;
441  if(gHarpoDebug>1)
442  Info("CompareHist","n = %d, W = %g, mean = %g, sigma = %g",n, W, mean, sigma);
443  // Double_t min = 0, max = 0, norm = 0, truncSum2 = 0;
444  // Double_t truncSum = HarpoTools::TruncSum(arr, 0.2, 0.8, min, max, norm, truncSum2,0.01);
445  delete arr;
446  // if(norm==0) return 0;
447 
448  // mean = truncSum/norm;
449  // sigma = truncSum2/norm - mean*mean;
450  // sigma = TMath::Sqrt(sigma);
451  // if(mean) sigma /= mean;
452  // if(gHarpoDebug>1)
453  // Info("CompareHist","** mean = %g, sigma = %g",mean, sigma);
454 
455 
456  return cov;
457 }
458 
459 
460 Double_t HarpoMatchingVertex::Covariance(TH1F* hX1, TH1F* hY1, TH1F* hX2, TH1F* hY2, TH1F* hoverlapX, TH1F* hoverlapY)
461 {
462 
463  Int_t nbins = hX1->GetXaxis()->GetNbins();
464  if(hY1->GetXaxis()->GetNbins() != nbins) return -1;
465 
466  Int_t n = 0;
467  Double_t W = 0, meanX = 0, sigX = 0, meanY = 0, sigY = 0, cov = 0;
468  Double_t qXold[2] = {-1000,-1000};
469  Double_t qYold[2] = {-1000,-1000};
470  for(Int_t bin = 1; bin<=nbins; bin++){
471  Double_t qX[2], qY[2];
472  qX[0] = hX1->GetBinContent(bin);
473  qY[0] = hY1->GetBinContent(bin);
474  qX[1] = hX2->GetBinContent(bin);
475  qY[1] = hY2->GetBinContent(bin);
476  Double_t dqX[2] = {-10000,-10000};
477  Double_t dqY[2] = {-10000,-10000};
478  for(Int_t i = 0; i<2; i++){
479  if(qXold[i]>0) dqX[i] = qX[i] - qXold[i];
480  if(qYold[i]>0) dqY[i] = qY[i] - qYold[i];
481  if(fPerfectMatching==3){
482  dqX[i] = qX[i];
483  dqY[i] = qY[i];
484  }
485  if(fPerfectMatching==4){
486  if(qXold[i]>0) dqX[i] = (qX[i] - qXold[i])/(qX[i]+qXold[i]);
487  if(qYold[i]>0) dqY[i] = (qY[i] - qYold[i])/(qY[i]+qYold[i]);
488  }
489  qXold[i] = qX[i];
490  qYold[i] = qY[i];
491  }
492  for(Int_t i = 0; i<2; i++){
493  if(qY[i]<=0) continue;
494  if(qX[i]<=0) continue;
495  if(dqX[i]==-10000) continue;
496  if(dqY[i]==-10000) continue;
497  Double_t w = 1;
498  if(hoverlapX->GetBinContent(bin)>0) w *= 0.1;
499  if(hoverlapY->GetBinContent(bin)>0) w *= 0.1;
500  // cov += qX[i]*qY[i];
501  // meanX += qX[i];
502  // sigX += qX[i]*qX[i];
503  // meanY += qY[i];
504  // sigY += qY[i]*qY[i];
505  cov += w*dqX[i]*dqY[i];
506  meanX += w*dqX[i];
507  sigX += w*dqX[i]*dqX[i];
508  meanY += w*dqY[i];
509  sigY += w*dqY[i]*dqY[i];
510  W += w;
511  n++;
512  }
513  }
514 
515  if(n==0) return 0;
516  meanX /= W;
517  sigX /= W;
518  sigX -= meanX*meanX;
519  meanY /= W;
520  sigY /= W;
521  sigY -= meanY*meanY;
522  cov /= W;
523  cov -= meanX*meanY;
524  cov /= TMath::Sqrt(sigX*sigY);
525 
526  return cov;
527 }
528 
529 
530 
531 void HarpoMatchingVertex::FillQvT(HarpoRecoVertex* vertex, TH1F* h1, TH1F* h2, TH1F* hoverlap, HarpoDccMap* m)
532 {
533 
534  Double_t xV = vertex->GetXvertex();
535  Double_t zV = vertex->GetZvertex();
536  Double_t pxV = vertex->GetPxVertex();
537  Double_t pzV = vertex->GetPzVertex();
538  Double_t thV = vertex->GetThetaVertex();
539  // Double_t cTh = TMath::Cos(thV*0.5);
540  // Double_t sTh = TMath::Sin(thV*0.5);
541  // Double_t px[2], pz[2];
542  // px[0] = pxV*cTh + pzV*sTh;
543  // pz[0] = pzV*cTh - pxV*sTh;
544  // px[1] = pxV*cTh - pzV*sTh;
545  // pz[1] = pzV*cTh + pxV*sTh;
546 
547  // Double_t fWidthX = 50, fWidthZ = 50;
548 
549  for(Int_t i=0;i<NADC;i++){ //Time bins
550  for(Int_t j=0;j<NCHAN;j++){ // Channels
551  Double_t q = m->GetData(j,i);
552  if(q<=0) continue;
553  // if(TMath::Abs(j-xV-px[0]/pz[0]*(i-zV))>fWidthX && TMath::Abs(j-xV-px[1]/pz[1]*(i-zV))>fWidthX ) continue;
554  // if(TMath::Abs(i-zV-pz[0]/px[0]*(j-xV))>fWidthZ && TMath::Abs(i-zV-pz[1]/px[1]*(j-xV))>fWidthZ) continue;
555  // Int_t track = 0;
556  // Double_t xi = xV + pxV/pzV*(i-zV);
557  // Double_t xi2 = xV + pxV/pzV*(i+1-zV);
558  // Double_t zj = zV + pzV/pxV*(j-xV);
559  // Double_t zj2 = zV + pzV/pxV*(j+1-xV);
560  Double_t r = 1; // r is the fraction of area in the sector containing the corner i,j
561  if(thV != 0 && (i-zV)*pxV-(j-xV)*pzV>0) r = 1. - r;
562  if(thV != 0 && TMath::Abs((i-zV)*pxV-(j-xV)*pzV)<1){
563  //Info("FillQvT","Overlap %g",(i-zV)*pxV-(j-xV)*pzV);
564  hoverlap->Fill(i,500);
565  }
566  h1->Fill(i,q*r);
567  h2->Fill(i,q*(1.-r));
568  }
569  }
570 
571  // for(Int_t i = 1; i<512; i++){
572  // Double_t q1 = h1->GetBinContent(i);
573  // Double_t q2 = h1->GetBinContent(i+1);
574  // h1->SetBinContent(i,q2-q1);
575  // q1 = h2->GetBinContent(i);
576  // q2 = h2->GetBinContent(i+1);
577  // h2->SetBinContent(i,q2-q1);
578  // }
579 
580 }
581 
582 
583 TH2F* HarpoMatchingVertex::Correl(TH1F* h1, TH1F* h2)
584 {
585 
586  Int_t nbins = h1->GetXaxis()->GetNbins();
587 
588  // TH2F* h = HistLogLog("h","",100,100,100000,100,100,100000);
589 
590  // for(Int_t i = 1; i<=nbins; i++)
591  // h->Fill(h1->GetBinContent(i),h2->GetBinContent(i));
592 
593  TH2F* h = new TH2F("h","",100,-10000,10000,100,-10000,10000);
594  Double_t q1old = -1000, q2old = -1000;
595  for(Int_t i = 1; i<=nbins; i++){
596  Double_t q1 = h1->GetBinContent(i);
597  Double_t q2 = h2->GetBinContent(i);
598  Double_t dq1 = -10000;
599  if(q1old>0) dq1 = q1-q1old;
600  Double_t dq2 = -10000;
601  if(q2old>0) dq2 = q2-q2old;
602  q1old = q1;
603  q2old = q2;
604  if(dq1==-10000) continue;
605  if(dq2==-10000) continue;
606  if(q1<0) continue;
607  if(q2<0) continue;
608  h->Fill(dq1,dq2);
609  }
610 
611  return h;
612 
613 
614 }
615 
616 
618 {
619 
620  fChooseNbins = 0;
621  fPerfectMatching = 3;
622  Long64_t perfectMatching = 0;
623  if ( !gHConfig->Lookup("matching.matchingMethod",perfectMatching) )
624  Info("Init","Use normal matching");
625  else{
626  fPerfectMatching = perfectMatching;
627  switch(fPerfectMatching){
628  case 0: Info("Init","Use normal matching"); break;
629  case 1: Info("Init","Use perfect matching (sim)"); break;
630  case 2: Info("Init","Use random matching"); break;
631  case 3: Info("Init","Use Q"); break;
632  case 4: Info("Init","Use (Q-Qold)/(Q+Qold)"); break;
633  default: Info("Init","Use normal matching"); break;
634  }
635  }
636 
637  // Initialise histograms here
638  for(Int_t i = 0; i<kMaxNvertex; i++){
639  for(Int_t j = 0; j<3; j++){
640  hQvT[i][j] = new TH1F(Form("hQvT_%d_%d",i,j),"",512,0,512);
641  }
642  hOverlap[i] = new TH1F(Form("hOverlap_%d",i),"",512,0,512);
643  }
644 
645  // hCov = HistLog("hCov","",100,0.001,100);
646  hCov = new TH1F("hCov","",100,0,10);
647  hCov2 = new TH2F("hCov2","",200,0,2,100,0,2);
648  hCov3 = new TH2F("hCov3","",200,-1,1,100,-1,1);
649 
650  fNtuple = new TNtupleD("fNtuple","","z0:x0:y0:px1:px2:py1:py2:cov:cov11:cov22:cov12:cov21:covSame:covSwitch:pxsim1:pysim1:pzsim1:pxsim2:pysim2:pzsim2:switchsim:switch:omegasim");
651  fNtuple->SetDirectory(0);
652 
653  }
654 
655 void HarpoMatchingVertex::Save(char * /* mode */)
656  {
657 
658  OpenHistFile("recoMatchingVertex");
659 
660  // Save histograms here
661  hCov->Write();
662  hCov2->Write();
663  hCov3->Write();
664  fNtuple->Write();
665  }
666 
667 
668 
669 void HarpoMatchingVertex::DisplayAnalysis(TRootEmbeddedCanvas* ecTab, TGListBox* infobox)
670 {
671  // in hrecomonitorgui: fills the canvas on the "Analysis" tab
672  TCanvas* c = ecTab->GetCanvas();
673  c->GetPad(1)->Delete();
674  c->GetPad(2)->Delete();
675  c->Divide(2);
676 
677  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
678  TClonesArray* vArray = reco->GetVertexArray();
679  Int_t nV = vArray->GetEntries();
680  Int_t nVp[2] = {0,0};
681  Int_t iV[2][nV];
682  for(Int_t i = 0; i<nV; i++){
683  HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(i);
684  Int_t plane = vertex->GetPlane();
685  iV[plane][nVp[plane]] = i;
686  nVp[plane]++;
687  }
688 
689 
690  // for(Int_t plane = 0; plane<2;plane++){
691  // TVirtualPad* cMap = c->GetPad(plane+1);
692  // cMap->Divide(1,nVp[plane]);
693  // for(Int_t i = 0; i<nVp[plane]; i++){
694  // hQvT[iV[plane][i]][1]->SetLineColor(kRed);
695  // hQvT[iV[plane][i]][2]->SetLineColor(kBlue);
696  // MakeNiceHisto(hQvT[iV[plane][i]][2],cMap->GetPad(i+1));
697  // hQvT[iV[plane][i]][1]->DrawCopy("same");
698  // }
699  // }
700 
701 
702 
703 
704  TLatex* latex = new TLatex();
705  latex->SetTextFont(132);
706  latex->SetTextAlign(12);
707  latex->SetTextSize(0.08);
708 
709 
710 
711  TClonesArray* v3dArray = reco->GetVertex3DArray();
712  Int_t nV3d = v3dArray->GetEntries();
713  c->GetPad(1)->Divide(4,nV3d+1);
714  c->GetPad(2)->Divide(4,nV3d+1);
715  MakeNiceHisto(hCov,c->GetPad(1)->GetPad(4*nV3d+1));
716  //c->GetPad(1)->GetPad(4*nV3d+1)->SetLogx();
717  MakeNiceHisto(hCov2,c->GetPad(2)->GetPad(4*nV3d+1),"col");
718  MakeNiceHisto(hCov3,c->GetPad(2)->GetPad(4*nV3d+2),"col");
719  // c->GetPad(2)->GetPad(2*nV3d+1)->SetLogx();
720  // c->GetPad(2)->GetPad(2*nV3d+1)->SetLogy();
721 
722  //infobox->RemoveAll();
723  infobox->NewEntry(Form("%d vertexes, %d 3D vertexes",nV,nV3d));
724  for(Int_t k = 0; k<nV3d; k++){
725  HarpoRecoVertex3D* vertex = (HarpoRecoVertex3D*)v3dArray->At(k);
726  infobox->NewEntry(Form("%d: (%g,%g,%g) 1:(%g,%g,%g) 2:(%g,%g,%g)",
727  k,
728  vertex->GetXvertex3D(),
729  vertex->GetYvertex3D(),
730  vertex->GetZvertex3D(),
731  vertex->GetPx1Vertex3D(),
732  vertex->GetPy1Vertex3D(),
733  vertex->GetPz1Vertex3D(),
734  vertex->GetPx2Vertex3D(),
735  vertex->GetPy2Vertex3D(),
736  vertex->GetPz2Vertex3D()));
737 
738  Int_t iVx = vertex->GetVertexX();
739  Int_t iVy = vertex->GetVertexY();
740  hQvT[iVy][1]->SetLineColor(kRed);
741  hQvT[iVy][2]->SetLineColor(kRed);
742  hQvT[iVx][1]->SetLineColor(kBlue);
743  hQvT[iVx][2]->SetLineColor(kBlue);
744  TH1F* hRatio11 = (TH1F*)hQvT[iVx][1]->Clone("hRatio11");
745  hRatio11->Divide(hQvT[iVy][1]);
746  TH1F* hRatio22 = (TH1F*)hQvT[iVx][2]->Clone("hRatio22");
747  hRatio22->Divide(hQvT[iVy][2]);
748  TH1F* hRatio12 = (TH1F*)hQvT[iVx][1]->Clone("hRatio12");
749  hRatio12->Divide(hQvT[iVy][2]);
750  TH1F* hRatio21 = (TH1F*)hQvT[iVx][2]->Clone("hRatio21");
751  hRatio21->Divide(hQvT[iVy][1]);
752 
753  // hRatio11->Scale(2000);
754  // hRatio22->Scale(2000);
755  // hRatio12->Scale(2000);
756  // hRatio21->Scale(2000);
757  // hRatio11->SetLineColor(kBlack);
758  // hRatio22->SetLineColor(kBlack);
759  // hRatio12->SetLineColor(kBlack);
760  // hRatio21->SetLineColor(kBlack);
761 
762  // // MakeNiceHisto(hQvT[iVx][1],c->GetPad(1)->GetPad(2*i+1));
763  // // hQvT[iVy][1]->DrawCopy("same");
764  // // //hRatio11->DrawCopy("same");
765  // // MakeNiceHisto(hQvT[iVx][2],c->GetPad(1)->GetPad(2*i+2));
766  // // hQvT[iVy][2]->DrawCopy("same");
767  // // //hRatio22->DrawCopy("same");
768  // // MakeNiceHisto(hQvT[iVx][1],c->GetPad(2)->GetPad(2*i+1));
769  // // hQvT[iVy][2]->DrawCopy("same");
770  // // //hRatio12->DrawCopy("same");
771  // // MakeNiceHisto(hQvT[iVx][2],c->GetPad(2)->GetPad(2*i+2));
772  // // hQvT[iVy][1]->DrawCopy("same");
773  // // //hRatio21->DrawCopy("same");
774 
775  TFile* outfile;
776  Bool_t save = kFALSE;
777  if(save) outfile = new TFile(Form("exampleMatching_%d.root",k),"recreate");
778 
779  Double_t mu, sig;
780  TH2F* hCorrel11 = Correl(hQvT[iVx][1],hQvT[iVy][1]);
781  if(save) hCorrel11->Write("hCorrel11");
782  // MakeNiceHisto(hCorrel11,c->GetPad(1)->GetPad(4*k+3),"col");
783  Double_t cov11 = CompareHist(hQvT[iVx][1],hQvT[iVy][1],hOverlap[iVx],hOverlap[iVy],mu,sig);
784  Double_t cov22 = CompareHist(hQvT[iVx][2],hQvT[iVy][2],hOverlap[iVx],hOverlap[iVy],mu,sig);
785  //cov11 = Covariance(hQvT[iVx][1],hQvT[iVy][1],hQvT[iVx][2],hQvT[iVy][2]);
786  TH2F* hCorrel22 = Correl(hQvT[iVx][2],hQvT[iVy][2]);
787  if(save) hCorrel22->Write("hCorrel22");
788  hCorrel22->Add(hCorrel11);
789  //latex->DrawLatex(500,5000,Form("%g",cov11));
790  TH2F* hCorrel12 = Correl(hQvT[iVx][1],hQvT[iVy][2]);
791  if(save) hCorrel12->Write("hCorrel12");
792  // MakeNiceHisto(hCorrel12,c->GetPad(2)->GetPad(4*k+3),"col");
793  Double_t cov12 = CompareHist(hQvT[iVx][1],hQvT[iVy][2],hOverlap[iVx],hOverlap[iVy],mu,sig);
794  Double_t cov21 = CompareHist(hQvT[iVx][2],hQvT[iVy][1],hOverlap[iVx],hOverlap[iVy],mu,sig);
795  //cov12 = Covariance(hQvT[iVx][1],hQvT[iVy][2],hQvT[iVx][2],hQvT[iVy][1]);
796  TH2F* hCorrel21 = Correl(hQvT[iVx][2],hQvT[iVy][1]);
797  if(save) hCorrel21->Write("hCorrel21");
798  hCorrel21->Add(hCorrel12);
799 
800  MakeNiceHisto(hCorrel22,c->GetPad(1)->GetPad(4*k+4),"col");
801  latex->DrawLatex(500,5000,Form("%g",TMath::Max(cov11,cov22)));
802  MakeNiceHisto(hCorrel21,c->GetPad(2)->GetPad(4*k+4),"col");
803  latex->DrawLatex(500,5000,Form("%g",TMath::Max(cov21,cov12)));
804 
805 
806  TVirtualPad* cMap[2];
807  cMap[0] = c->GetPad(1);
808  cMap[1] = c->GetPad(2);
809  for(Int_t plane =0; plane<2; plane++){
810  HarpoDccMap* m = fEvt->GetDccMap(plane);
811  HarpoRecoVertex* vertex;
812  if(plane == XPLANE) vertex = (HarpoRecoVertex*)vArray->At(iVx);
813  if(plane == YPLANE) vertex = (HarpoRecoVertex*)vArray->At(iVy);
814  Double_t xV = vertex->GetXvertex();
815  Double_t zV = vertex->GetZvertex();
816  Double_t pxV = vertex->GetPxVertex();
817  Double_t pzV = vertex->GetPzVertex();
818  Double_t thV = vertex->GetThetaVertex();
819  // Double_t cTh = TMath::Cos(thV*0.5);
820  // Double_t sTh = TMath::Sin(thV*0.5);
821  // Double_t px[2], pz[2];
822  // px[0] = pxV*cTh + pzV*sTh;
823  // pz[0] = pzV*cTh - pxV*sTh;
824  // px[1] = pxV*cTh - pzV*sTh;
825  // pz[1] = pzV*cTh + pxV*sTh;
826  TH2F* h1 = new TH2F(Form("h1_%d",plane),"",512,0,512,288,0,288);
827  TH2F* h2 = new TH2F(Form("h2_%d",plane),"",512,0,512,288,0,288);
828  // TH1F* p1 = new TH1F("p1","",512,0,512);
829  // TH1F* p2 = new TH1F("p2","",512,0,512);
830  for(Int_t i=0;i<NADC;i++){ //Time bins
831  for(Int_t j=0;j<NCHAN;j++){ // Channels
832  Double_t q = m->GetData(j,i);
833  if(q<=0) continue;
834  // Exclude data not belonging to tracks
835  // if(TMath::Abs(j-xV-px[0]/pz[0]*(i-zV))>5 &&
836  // TMath::Abs(j-xV-px[1]/pz[1]*(i-zV))>5 &&
837  // TMath::Abs(i-zV-pz[0]/px[0]*(j-xV))>10 &&
838  // TMath::Abs(i-zV-pz[1]/px[1]*(j-xV))>10) continue;
839  // Int_t track = 0;
840  // Double_t xi = xV + pxV/pzV*(i-zV);
841  // Double_t xi2 = xV + pxV/pzV*(i+1-zV);
842  // Double_t zj = zV + pzV/pxV*(j-xV);
843  // Double_t zj2 = zV + pzV/pxV*(j+1-xV);
844  Double_t r = 1; // r is the fraction of area in the sector containing the corner i,j
845  if(thV != 0 && (i-zV)*pxV-(j-xV)*pzV>0) r = 1. - r;
846  h1->SetBinContent(i+1,j+1,q*r);
847  h2->SetBinContent(i+1,j+1,q*(1.-r));
848  // p1->Fill(i,q*r);
849  // p2->Fill(i,q*(1.-r));
850  }
851  }
852  h1->SetLineColor(kGreen);
853  h2->SetLineColor(kMagenta);
854  MakeNiceHisto(h1,cMap[plane]->GetPad(4*k+1),"box");
855  h2->DrawCopy("box same");
856  latex->DrawLatex(50,200,Form("%d",iV[plane][k]));
857 
858  // if(plane==0){
859  // p1->SetLineColor(kRed);
860  // MakeNiceHisto(p1,cMap[0]->GetPad(4*k+2));
861  // p2->SetLineColor(kRed);
862  // MakeNiceHisto(p2,cMap[0]->GetPad(4*k+3));
863  // p1->SetLineColor(kRed);
864  // MakeNiceHisto(p1,cMap[1]->GetPad(4*k+2));
865  // p2->SetLineColor(kRed);
866  // MakeNiceHisto(p2,cMap[1]->GetPad(4*k+3));
867  // }else{
868  // cMap[0]->cd(4*k+2);
869  // p1->SetLineColor(kBlue);
870  // p1->DrawCopy("same");
871  // cMap[0]->cd(4*k+3);
872  // p2->SetLineColor(kBlue);
873  // p2->DrawCopy("same");
874  // cMap[1]->cd(4*k+2);
875  // p2->SetLineColor(kBlue);
876  // p2->DrawCopy("same");
877  // cMap[1]->cd(4*k+3);
878  // p1->SetLineColor(kBlue);
879  // p1->DrawCopy("same");
880  // }
881 
882  if(save) h1->Write();
883  if(save) h2->Write();
884 
885  h1->Delete();
886  h2->Delete();
887  // p1->Delete();
888  // p2->Delete();
889 
890  }
891  if(save) hQvT[iVx][1]->Write("hQvT_X_1");
892  if(save) hQvT[iVx][2]->Write("hQvT_X_2");
893  if(save) hQvT[iVy][1]->Write("hQvT_Y_1");
894  if(save) hQvT[iVy][2]->Write("hQvT_Y_2");
895 
896  if(save) outfile->Close();
897 
898  hQvT[iVx][1]->SetLineColor(kRed);
899  hQvT[iVx][2]->SetLineColor(kRed);
900  hQvT[iVy][1]->SetLineColor(kBlue);
901  hQvT[iVy][2]->SetLineColor(kBlue);
902 
903  hOverlap[iVx]->SetFillColor(kRed);
904  hOverlap[iVy]->SetFillColor(kBlue);
905  hOverlap[iVx]->SetFillStyle(3004);
906  hOverlap[iVy]->SetFillStyle(3005);
907 
908  MakeNiceHisto(hQvT[iVx][1],cMap[0]->GetPad(4*k+2));
909  hQvT[iVy][1]->DrawCopy("same");
910  hOverlap[iVx]->DrawCopy("same");
911  hOverlap[iVy]->DrawCopy("same");
912 
913  MakeNiceHisto(hQvT[iVx][2],cMap[0]->GetPad(4*k+3));
914  hQvT[iVy][2]->DrawCopy("same");
915  hOverlap[iVx]->DrawCopy("same");
916  hOverlap[iVy]->DrawCopy("same");
917 
918  MakeNiceHisto(hQvT[iVx][1],cMap[1]->GetPad(4*k+2));
919  hQvT[iVy][2]->DrawCopy("same");
920  hOverlap[iVx]->DrawCopy("same");
921  hOverlap[iVy]->DrawCopy("same");
922 
923  MakeNiceHisto(hQvT[iVx][2],cMap[1]->GetPad(4*k+3));
924  hQvT[iVy][1]->DrawCopy("same");
925  hOverlap[iVx]->DrawCopy("same");
926  hOverlap[iVy]->DrawCopy("same");
927  }
928 
929 
930 
931 
932  // for(Int_t plane = 0; plane<2;plane++){
933  // TVirtualPad* cMap = c->GetPad(plane+1);
934  // HarpoDccMap* m = fEvt->GetDccMap(plane);
935  // // cMap->Divide(2,nVp[plane]);
936  // for(Int_t k = 0; k<nVp[plane]; k++){
937  // HarpoRecoVertex* vertex = (HarpoRecoVertex*)vArray->At(iV[plane][k]);
938  // Double_t xV = vertex->GetXvertex();
939  // Double_t zV = vertex->GetZvertex();
940  // Double_t pxV = vertex->GetPxVertex();
941  // Double_t pzV = vertex->GetPzVertex();
942  // Double_t thV = vertex->GetThetaVertex();
943  // Double_t cTh = TMath::Cos(thV*0.5);
944  // Double_t sTh = TMath::Sin(thV*0.5);
945  // Double_t px[2], pz[2];
946  // px[0] = pxV*cTh + pzV*sTh;
947  // pz[0] = pzV*cTh - pxV*sTh;
948  // px[1] = pxV*cTh - pzV*sTh;
949  // pz[1] = pzV*cTh + pxV*sTh;
950  // TH2F* h1 = new TH2F("h1","",512,0,512,288,0,288);
951  // TH2F* h2 = new TH2F("h2","",512,0,512,288,0,288);
952  // TH1F* p1 = new TH1F("p1","",512,0,512);
953  // TH1F* p2 = new TH1F("p2","",512,0,512);
954  // for(Int_t i=0;i<NADC;i++){ //Time bins
955  // for(Int_t j=0;j<NCHAN;j++){ // Channels
956  // Double_t q = m->GetData(j,i);
957  // if(q<=0) continue;
958  // // Exclude data not belonging to tracks
959  // if(TMath::Abs(j-xV-px[0]/pz[0]*(i-zV))>5 &&
960  // TMath::Abs(j-xV-px[1]/pz[1]*(i-zV))>5 &&
961  // TMath::Abs(i-zV-pz[0]/px[0]*(j-xV))>10 &&
962  // TMath::Abs(i-zV-pz[1]/px[1]*(j-xV))>10) continue;
963  // Int_t track = 0;
964  // Double_t xi = xV + pxV/pzV*(i-zV);
965  // Double_t xi2 = xV + pxV/pzV*(i+1-zV);
966  // Double_t zj = zV + pzV/pxV*(j-xV);
967  // Double_t zj2 = zV + pzV/pxV*(j+1-xV);
968  // Double_t r = 1; // r is the fraction of area in the sector containing the corner i,j
969  // if(thV != 0 && (i-zV)*pxV-(j-xV)*pzV>0) r = 1. - r;
970  // h1->SetBinContent(i+1,j+1,q*r);
971  // h2->SetBinContent(i+1,j+1,q*(1.-r));
972  // p1->Fill(i,q*r);
973  // p2->Fill(i,q*(1.-r));
974  // }
975  // }
976  // h1->SetLineColor(kRed);
977  // h2->SetLineColor(kBlue);
978  // MakeNiceHisto(h1,cMap->GetPad(4*k+1),"box");
979  // h2->DrawCopy("box same");
980  // latex->DrawLatex(50,200,Form("%d",iV[plane][k]));
981 
982 
983  // p1->SetLineColor(kRed);
984  // p2->SetLineColor(kBlue);
985  // MakeNiceHisto(p1,cMap->GetPad(4*k+2));
986  // p2->DrawCopy("same");
987 
988  // h1->Delete();
989  // h2->Delete();
990  // p1->Delete();
991  // p2->Delete();
992  // }
993 
994  // }
995 
996 
997 
998 
999  c->Update();
1000 }
1001 
1002 
1003 void HarpoMatchingVertex::ConfigFrame(TGMainFrame* fMain, Int_t id)
1004 {
1005  // in hrecomonitorgui: Creates a popup window for analysis configuration
1006  // You must define SetConfig() properly
1007 
1008  UInt_t xsize = 200;
1009  UInt_t ysize2 = 20;
1010  UInt_t ysize = 10*ysize2;
1011  TGTransientFrame* main = new TGTransientFrame(gClient->GetRoot(), fMain, xsize, ysize);
1012  main->Connect("CloseWindow()", "HarpoRecoMonitorGui", main, "CloseWindow()");
1013  main->DontCallClose(); // to avoid double deletions.
1014 
1015  // use hierarchical cleaning
1016  main->SetCleanup(kDeepCleanup);
1017 
1018  TGVerticalFrame* frame = new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);
1019 
1020  // Object layout options
1021  TGLayoutHints* fLayout1 = new TGLayoutHints(kLHintsLeft ,5,5,5,5);
1022  TGLayoutHints* fLayout2 = new TGLayoutHints(kLHintsRight ,5,5,5,5);
1023  TGLayoutHints* fLayout3 = new TGLayoutHints(kLHintsTop | kLHintsExpandX ,5,5,5,5);
1024  TGLayoutHints* fLayout4 = new TGLayoutHints(kLHintsExpandX | kLHintsExpandY,5,5,5,5);
1025 
1026  //________ DO NOT MODIFY ABOVE _____________________________
1027 
1028 
1029 
1030 
1031  // Title of the analysis
1032  TGLabel* fAnalysisLabel = new TGLabel(frame, "Template analysis");
1033 
1034  // Create a line for choosing value of parameter
1035  TGHorizontalFrame* nBinsFrame = new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
1036  TGLabel* nBinsLabel = new TGLabel(nBinsFrame, "fNbins");
1037  fChooseNbins = new TGNumberEntry(nBinsFrame);
1038  fChooseNbins->SetNumber(0);
1039 
1040  main->AddFrame(frame,fLayout4);
1041  frame->AddFrame(fAnalysisLabel,fLayout3);
1042  frame->AddFrame(nBinsFrame,fLayout3);
1043  nBinsFrame->AddFrame(nBinsLabel,fLayout1);
1044  nBinsFrame->AddFrame(fChooseNbins,fLayout2);
1045 
1046 
1047 
1048  //________ DO NOT MODIFY BELOW _____________________________
1049  // Button to validate configuration
1050  TGTextButton* setConf = new TGTextButton(frame,"Save Config",id);
1051  setConf->Associate(fMain);
1052 
1053  frame->AddFrame(setConf,fLayout3);
1054 
1055  main->MapSubwindows();
1056  main->MapWindow();
1057  main->Resize();
1058  return;
1059 }
1060 
1061 
1062 
1064 {
1065  // Update the configuration according to the values in the popup window
1066 
1067  if(!fChooseNbins) return;
1068 
1069  // fNbins = fChooseNbins->GetNumber();
1070 }
Double_t GetZvertex3D()
Double_t GetYvertex3D()
void FillQvT(HarpoRecoVertex *vertex, TH1F *h1, TH1F *h2, TH1F *hoverlap, HarpoDccMap *m)
Redefine empty default.
Double_t CompareHist(TH1F *hX, TH1F *hY, TH1F *hoverlapX, TH1F *hoverlapY, Double_t &mean, Double_t &sigma)
void AddVertex3D(HarpoRecoVertex3D *val)
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
Double_t Covariance(TH1F *hX1, TH1F *hY1, TH1F *hX2, TH1F *hY2, TH1F *hoverlapX, TH1F *hoverlapY)
#define XPLANE
Definition: HarpoConfig.h:24
void DisplayAnalysis(TRootEmbeddedCanvas *ecTab, TGListBox *infobox)
#define NCHAN
Definition: HarpoDccMap.h:16
void Save(char *mode=NULL)
Double_t GetPz1Vertex3D()
Double_t GetPy2Vertex3D()
void print()
Overloaded method which do all job.
Double_t GetZvertex()
Double_t GetPz()
Definition: TpcSimMCEvent.h:32
HarpoMatchingVertex.
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
int main(int argc, char **argv)
Definition: drawAnalyse.cxx:25
Double_t GetThetaVertex()
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
Definition: HarpoDetSet.h:33
TGNumberEntry * fChooseNbins
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
2D vertex object, containing position, angle and associated track numbers, and quality info ...
TH2F * Correl(TH1F *h1, TH1F *h2)
void ConfigFrame(TGMainFrame *fMain, Int_t id)
virtual void print()
Double_t GetPy()
Definition: TpcSimMCEvent.h:31
FullEvent Header not scecific to the detectors The class is ....
Double_t GetPzVertex()
TClonesArray * GetVertexArray()
A class store HARPO row DCC event data and header. End provide access metods to the row data...
Definition: HarpoSimEvent.h:24
Double_t GetPx1Vertex3D()
virtual void print()
TFile * OpenHistFile(const char *ananame)
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
HarpoDetEvent * GetDetEvent(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:93
TH1F * hOverlap[kMaxNvertex]
Double_t GetPy1Vertex3D()
Data from Keller temperuture and pressure sensors.
Definition: HarpoDet.h:22
TH1F * hQvT[kMaxNvertex][3]
void ResetVertex3DArray()
TVirtualPad * MakeNiceHisto(TH1 *hist, TVirtualPad *c1, const char *opt, Bool_t copy)
Int_t GetNtracks()
TpcSimMCTrack * GetTrack(Int_t i)
Definition: TpcSimMCEvent.h:91
#define NADC
Definition: HarpoDccMap.h:18
Long64_t gHarpoDebug
Definition: HarpoDebug.cxx:9
#define YPLANE
Definition: HarpoConfig.h:25
TpcSimMCEvent * GetMCEvent()
Definition: HarpoSimEvent.h:57
Bool_t Lookup(const char *path, Bool_t &val)
Lookup function for scalar values.
TClonesArray * GetVertex3DArray()
static const Int_t kMaxNvertex
HarpoEventHeader * GetHeader()
Definition: HarpoEvent.cxx:80
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
ULong_t nEvents
Definition: HarpoAnalyse.h:75
Double_t GetPx2Vertex3D()
Double_t GetPz2Vertex3D()
Double_t GetPxVertex()
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
HarpoConfig * gHConfig
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
Double_t GetPx()
Definition: TpcSimMCEvent.h:30
Double_t GetXvertex()
R__EXTERN HarpoDetSet * gHDetSet
Definition: HarpoDetSet.h:71
Double_t GetXvertex3D()
3D vertex object, containing position, angle and associated 2D vertexes, and quality info ...