HARPO  5.1.1
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
HarpoKalmanTracking.cxx
Go to the documentation of this file.
1 //
2 // File HarpoKalmanTracking.cxx
3 //
11 #include "HarpoKalmanTracking.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 
20 #include "TFile.h"
21 #include "TStyle.h"
22 #include "TCanvas.h"
23 #include "TLatex.h"
24 #include "TGraph.h"
25 #include "TF1.h"
26 #include "TArrow.h"
27 #include "TLinearFitter.h"
28 #include "TMath.h"
29 #include "TSystem.h"
30 
31 #include <cstdlib>
32 #include <cstring>
33 #include <cassert>
34 #include <fstream>
35 #include <iostream>
36 
37 
38 ClassImp(HarpoKalmanTracking);
39 
41  {
42  unsigned int nd; // number of detectors
44 
45  assert(hdr != NULL);
46  hdr->print();
47 
48  for (nd = 0; nd < gkNDetectors; nd++) {
49  // if (fEvt->isdataExist(nd)) {
50  HarpoDccMap *plane = fEvt->GetDccMap(nd);
51  if (plane != NULL )plane->print();
52  }
53  }
54 
56 {
57  // Bool_t drawEvent = kFALSE;
58  nEvents++;
59  //if(gHarpoDebug>0)
60  Info("process","Event %ld",nEvents);
61 
62 
63 
64 
65 
66  for(Int_t i = 0; i< NTRACK; i++)
67  fNclTrack[i] = 0;
68  if(!fEvt->GetRecoEvent()) // Reconstructed data (clusters, tracks, matched tracks)
69  return;
70  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
71  reco->SetTrackType(KALMAN);
72  // Clusters in the event (both X and Y directions)
73  TClonesArray* clArray = reco->GetClustersArray();
74  // Int_t nCl = clArray->GetEntries();
75 
76 
77  // TMatrixD Xorig(4,1);
78  Double_t inc = fResKalman; // Space resolution
79  Double_t theta = fScatKalman; // Multiple scattering
80  TMatrixD Corig(4,4);
81  for(int a = 0;a<4;a++){
82  for(int b =0;b<4;b++){
83  if(a==b){
84  if(a<2)
85  Corig[a][b]=inc*inc; //initialisation de C
86  else
87  Corig[a][b]=(2*inc)*(2*inc)+theta*theta;
88  }
89  else
90  Corig[a][b]=0;
91  }
92  }
93 
94  for(Int_t plane = 0; plane<2; plane++){
95  if(!InitPlane(clArray,plane)) continue;
96 
97  Int_t counter = 0, color = 51, nClLeftOld = 0, nTry = 0;
98  while(1){
99  if(counter>50){
100  //Warning("process","Too many loops on plane %d",plane);
101  break;
102  }
103  if(gHarpoDebug>1)
104  Info("process","Look for seeds on plane %d %d",plane,counter);
105  counter++;
106 
107  Int_t iMin = NCHAN, iMax = 0, jMin = NADC, jMax = 0;
108  Int_t nClLeft = GetMapEdges(clArray,plane,iMin,iMax,jMin,jMax);
109  if(nClLeft<10){
110  if(gHarpoDebug>0)
111  Info("process","Only %d clusters left",nClLeft);
112  break;
113  }
114  if(nClLeft==nClLeftOld){
115  if(gHarpoDebug>0)
116  Info("process","No new tracks, %d clusters left",nClLeft);
117  break;
118  }
119  nClLeftOld = nClLeft;
120 
121 
122 
123 
124  for(Int_t side = 0; side<4; side++){ // Loop over starting map side
125  Int_t type = -1, index0 = 0, index1 = 0, imax0 = 0, imax1 = 0, sign = 0;
126  switch(side){
127  case 0: type = TCLUSTER; index0 = jMax; sign = -1; break;
128  case 1: type = CCLUSTER; index0 = iMax; sign = -1; break;
129  case 2: type = CCLUSTER; index0 = iMin; sign = +1; break;
130  case 3: type = TCLUSTER; index0 = jMin; sign = +1; break;
131  }
132  // if(side == 1 || side == 2) gHarpoDebug = 3;
133  // else continue;//gHarpoDebug = 0;
134  index1 = index0 + sign;
135  if(type == TCLUSTER){
136  if(index0<0 || index0>=NADC) continue;
137  while(NTcl[index1]<=0 && TMath::Abs(index1-index0) < 50) index1 += sign;
138  if(index1<0 || index1>=NADC) continue;
139  imax0 = NTcl[index0];
140  imax1 = NTcl[index1];
141  }else{
142  if(index0<0 || index0>=NCHAN) continue;
143  while(NCcl[index1]<=0 && TMath::Abs(index1-index0) < 50) index1 += sign;
144  if(index1<0 || index1>=NCHAN) continue;
145  imax0 = NCcl[index0];
146  imax1 = NCcl[index1];
147  }
148 
149 
150  if(TMath::Abs(index1-index0) >= 10){
151  for(Int_t i0 = 0; i0<imax0; i0++){
152  if(i0 == 20) break;
153  Int_t icl0;
154  if(type == CCLUSTER) icl0 = Ccl[index0][i0];
155  else icl0 = Tcl[index0][i0];
156  // used[icl0] = 3;
157  used[icl0] = -2000;
158  }
159  if(gHarpoDebug>1) Info("process","side %d: i0 = %d, i1 = %d",side,index0,index1);
160  continue;
161  }
162 
163 
164  if(gHarpoDebug>1)
165  Info("process","side %d => index = %d %d (%d %d)",side,index0,index1,imax0,imax1);
166  for(Int_t i0 = 0; i0<imax0; i0++){ // Loop over starting cluster pairs
167  if(i0 == 20) break;
168  Int_t icl0;
169  if(type == CCLUSTER) icl0 = Ccl[index0][i0];
170  else icl0 = Tcl[index0][i0];
171  if(gHarpoDebug>1) Info("process","used[%d] = %d",icl0,used[icl0]);
172  // if(used[icl0]) continue;
173  // used[icl0] = 3;
174  if(used[icl0] != -1000) continue;
175  used[icl0] = -2000;
176  for(Int_t i1 = 0; i1<imax1; i1++){
177  Int_t icl1;
178  if(type == CCLUSTER) icl1 = Ccl[index1][i1];
179  else icl1 = Tcl[index1][i1];
180  if(gHarpoDebug>1) Info("process","used[%d] = %d",icl1,used[icl1]);
181  // if(used[icl1]) continue;
182  if(used[icl1] != -1000) continue;
183  //used[icl1] = 3;
184 
185  //if(nTry==500){
186  // Warning("process","TOO MANY TRACK SEEDS");
187  // return;
188  //}
189  if(gHarpoDebug>1) Info("process","nTry = %d",nTry);
190  FindTrack(clArray,icl0,icl1,Corig,plane,color);
191  nTry++;
192 
193  }
194  }
195  }
196 
197  }
198  SpliceTracks(plane);
199  //SpliceTracksNew(plane);
200 
201  if(fHist[plane]){
202  fHist[plane]->GetXaxis()->SetRangeUser(0,511);
203  fHist[plane]->GetYaxis()->SetRangeUser(0,304);
204  fCanvas[plane]->Update();
205  }
206  }
207 }
208 
209 Int_t HarpoKalmanTracking::InitPlane(TClonesArray* clArray, Int_t plane)
210 {
211  HarpoKalman::InitPlane(clArray, plane);
212  for(Int_t itr = 0; itr<NTRACK; itr++){
213  fId[itr][0] = -1;
214  fId[itr][1] = -1;
215  }
216  return 1;
217 }
218 
219 
220 
221 Int_t HarpoKalmanTracking::GetMapEdges(TClonesArray* clArray, Int_t plane, Int_t &iMin, Int_t &iMax, Int_t &jMin, Int_t &jMax)
222 {
223 
224  Int_t nClLeft = 0;
225  iMin = NCHAN+1; iMax = -1; jMin = NADC+1; jMax = -1;
226  Int_t nCl = clArray->GetEntries();
227  for(Int_t icl = 0; icl<nCl; icl++){
228  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
229  // if(cluster->GetIdClusterTrack(0)!=-1) continue;
230  // if(used[icl]) continue;
231  if(used[icl] != -1000) continue;
232  Int_t index = cluster->GetIndex();
233  if(cluster->GetQuality() != 0) continue;
234  if(cluster->GetPlane() != plane) continue;
235  if(cluster->GetQ() < fQmin) continue;
236  if(cluster->GetType() == TCLUSTER && cluster->GetWidth() < fWidthMinT) continue;
237  if(cluster->GetType() == CCLUSTER && cluster->GetWidth() < fWidthMinC) continue;
238  if(cluster->GetWidth() > fWidthMax) continue;
239  if(cluster->GetType() == TCLUSTER){
240  if(index>jMax) jMax = index;
241  if(index<jMin) jMin = index;
242  }
243  if(cluster->GetType() == CCLUSTER){
244  if(index>iMax) iMax = index;
245  if(index<iMin) iMin = index;
246  }
247  nClLeft++;
248  }
249 
250  if(gHarpoDebug>1)
251  Info("GetMapEdges","%d %d %d %d (%d)",iMin, iMax, jMin, jMax, nClLeft);
252 
253  return nClLeft;
254 }
255 
256 void HarpoKalmanTracking::FindTrack(TClonesArray* clArray, Int_t icl0, Int_t icl1, TMatrixD Corig, Int_t plane, Int_t &color)
257 {
258 
259  if(gHarpoDebug>2)
260  Info("FindTrack","%d %d",icl0,icl1);
261  if(fNtr[plane]>=NTRACK) return;
262 
263 
264  if(fGraphs[plane][fNtr[plane]] != 0){
265  fGraphs[plane][fNtr[plane]]->Delete();
266  fGraphs[plane][fNtr[plane]] = new TGraph();
267  }
268 
269  HarpoRecoClusters* cluster0 = (HarpoRecoClusters*)clArray->At(icl0);
270  HarpoRecoClusters* cluster1 = (HarpoRecoClusters*)clArray->At(icl1);
271  Int_t type = cluster0->GetType();
272  if(type != 0 && type != 1) return;
273  // if(cluster0->GetIdClusterTrack(0)!=-1 && cluster1->GetIdClusterTrack(0)!=-1) continue;
274  Double_t x0 = cluster0->GetMean();
275  Double_t x1 = cluster1->GetMean();
276  Double_t index0 = cluster0->GetIndex();
277  Double_t index1 = cluster1->GetIndex();
278  Int_t q0 = cluster0->GetQ();
279  // Int_t itr = fNtr[plane];
280  if(TMath::Abs(x0-x1)>15) return;
281  TMatrixD Xorig(4,1);
282  Xorig[type][0] = x0;
283  Xorig[1-type][0] = index0;
284  Xorig[2+type][0] = x1-x0;
285  Xorig[3-type][0] = index1-index0;
286  if(gHarpoDebug>2)
287  Info("process","FindNext(%f,%f,%f,%f,%d)",Xorig[0][0],Xorig[1][0],Xorig[2][0],Xorig[3][0],fNtr[plane]);
288  if(TMath::Abs(Xorig[2+type][0]/Xorig[3-type][0])>1.5) return; // Angular cut for each cluster type
289 
290  if(fCanvas[plane]){
291  TArrow* arrow = new TArrow();
292  arrow->SetLineColor(kBlack);
293  arrow->SetFillColor(color);
294  fCanvas[plane]->cd();
295  arrow->DrawArrow(Xorig[0][0],Xorig[1][0],Xorig[0][0]+Xorig[2][0],Xorig[1][0]+Xorig[3][0],0.01);
296  fHist[plane]->GetXaxis()->SetRangeUser(Xorig[0][0] - 10,Xorig[0][0] + 10);
297  fHist[plane]->GetYaxis()->SetRangeUser(Xorig[1][0] - 10,Xorig[1][0] + 10);
298  fCanvas[plane]->Update();
299  }
300 
301 
302  // TArrayI* arr = new TArrayI(1000);
303  // arr->AddAt(icl0+1,0); // Keep 0 and 1 for start and end
304  // arr =
305  // FindNext(Xorig,Corig,fNtr[plane],clArray,arr, plane,color,3001,0);
306  // FindNext(Xorig,Corig,fNtr[plane],clArray,0, plane,color,3001,0);
307  Double_t angleorig = TMath::ATan2(Xorig[2][0],Xorig[3][0]);
308  FindNext(Xorig,Corig,angleorig,fNtr[plane],clArray,0,1,plane,color,3001,0,q0);
309  if(gHarpoDebug>1)
310  Info("process","Found new track %d",fNtr[plane]);
311  // Int_t nClTr = AddTrack(arr,clArray,plane);
312  Int_t nClTr = AddTrack(clArray,plane);
313  if(nClTr>10)
314  color += 4;
315  // delete arr;
316  if(gHarpoDebug>1)
317  Info("process","Track %d: %d clusters",fNtr[plane]-1,nClTr);
318  // fNtr[plane]++;
319 
320 
321 }
322 
323 
324 
325 
326 //Int_t HarpoKalmanTracking::AddTrack(Int_t iTr, TClonesArray* clArray, Int_t plane)
327 //Int_t HarpoKalmanTracking::AddTrack(TArrayI* array, TClonesArray* clArray, Int_t plane)
328 Int_t HarpoKalmanTracking::AddTrack(TClonesArray* clArray, Int_t plane)
329 {
330 
331  Int_t n = 0;
332  if(fNtr[plane]<0 || fNtr[plane]>=NTRACK){
333  Warning("AddTrack","Wrong track number %d",fNtr[plane]);
334  return -10000;
335  }
336 
337  TMatrixD X(4,1);
338  Int_t index = (fStartIndex[fNtr[plane]]+9)%10;
339  //Info("AddTrack","%d %d %d",plane,fNtr[plane],index);
340  X[0][0] = fStartPointZ[fNtr[plane]][index] + fStartDirZ[fNtr[plane]][index];
341  X[1][0] = fStartPointX[fNtr[plane]][index] + fStartDirX[fNtr[plane]][index];
342  X[2][0] = -fStartDirZ[fNtr[plane]][index];
343  X[3][0] = -fStartDirX[fNtr[plane]][index];
344 
345  if(gHarpoDebug>1)
346  Info("AddTrack","New seed: (%g, %g, %g, %g)",X[0][0],X[1][0],X[2][0],X[3][0]);
347 
348  TMatrixD Corig(4,4);
349  for(int a = 0;a<4;a++){
350  for(int b =0;b<4;b++){
351  if(a==b){
352  if(a<2)
353  Corig[a][b]=fResKalman*fResKalman; //initialisation de C
354  else
355  Corig[a][b]=(2*fResKalman)*(2*fResKalman)+fScatKalman*fScatKalman;
356  }
357  else
358  Corig[a][b]=0;
359  }
360  }
361  TArrayI* arr = new TArrayI(2000);
362  Double_t angleorig = TMath::ATan2(X[2][0],X[3][0]);
363  arr = FindNextClosest(X, Corig, angleorig,fNtr[plane], clArray, arr, 1,plane, kGray, 3002, 2, -1);
364 
365  for(Int_t i = 0; i<2000; i++){
366  if(arr->At(i)-1<0) break;
367  n++;
368  }
369  //Info("AddTrack","Ncl = %d",n);
370  if(n<fNclMin){
371  delete arr;
372  return -n;
373  }
374 
375  if(fGraphs[plane][fNtr[plane]])
376  fGraphs[plane][fNtr[plane]]->Delete();
377  fGraphs[plane][fNtr[plane]] = new TGraph();
378  if(plane == XPLANE)
379  fGraphs[plane][fNtr[plane]]->SetName(Form("TrackX%d",fNtr[plane]));
380  else
381  fGraphs[plane][fNtr[plane]]->SetName(Form("TrackY%d",fNtr[plane]));
382 
383  for(Int_t i = 0; i<n; i++){
384  Int_t icl = arr->At(i) - 1;
385  //Info("AddTarack","icl = %d",icl);
386  if(icl<0) break;
387  if(icl>=4000) continue;
388  fTrack[icl] = fNtr[plane];
389  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
390  Double_t mean = cluster->GetMean();
391  Int_t index = cluster->GetIndex();
392  if(cluster->GetQuality() != 0) continue;
393  //cluster->AddIdClusterTrack(fNtr[plane]);
394  if(cluster->GetType() == TCLUSTER){
395  fGraphs[plane][fNtr[plane]]->SetPoint(fGraphs[plane][fNtr[plane]]->GetN(),index,mean);
396 
397  }
398  if(cluster->GetType() == CCLUSTER){
399  fGraphs[plane][fNtr[plane]]->SetPoint(fGraphs[plane][fNtr[plane]]->GetN(),mean,index);
400  }
401  }
402 
403  delete arr;
404 
405 
406  // if(gHarpoDebug>1)
407  // Info("AddTrack","Added track %d (plane %d) (%d cl) start(%g, %g) end(%g, %g)",fNtr[plane],plane,n,fStartPointX[fNtr[plane]][fStartIndex[fNtr[plane]]%10],fStartPointZ[fNtr[plane]][fStartIndex[fNtr[plane]]%10],fEndPointX[fNtr[plane]][fEndIndex[fNtr[plane]]%10],fEndPointZ[fNtr[plane]][fEndIndex[fNtr[plane]]%10]);
408 
409  fNtr[plane]++;
410 
411  return n;
412 }
413 
414 
415 
417 {
418  // return;
419 
420  // Double_t fDmin = 5, fThetaMin = -0.5;
421  Double_t fThetaMin = -0.5;
422  Double_t fDminForward = 20, fDminTransverse = 2;
423  // Double_t fDmin2 = fDmin*fDmin;
424  Int_t nSplice = 0;
425  Int_t fOffset = 0;
426  if(gHarpoDebug>1)
427  Info("SpliceTracks","plane %d: Ntracks = %d",plane,fNtr[plane]);
428 
429  Double_t x[fNtr[plane]][2], z[fNtr[plane]][2], th[fNtr[plane]][2];
430  Double_t px[fNtr[plane]][2], pz[fNtr[plane]][2];
431  Int_t spliced[fNtr[plane]];
432  GetQtracks(plane);
433  for(Int_t itr = 0; itr<fNtr[plane]; itr++){
434  spliced[itr] = itr;
435  if(fGraphs[plane][itr]->GetN()<=fOffset+1) continue;
436  fGraphs[plane][itr]->GetPoint(0 + fOffset,z[itr][0],x[itr][0]);
437  Double_t xTmp, zTmp;
438  fGraphs[plane][itr]->GetPoint(1 + fOffset,xTmp,zTmp);
439  px[itr][0] = xTmp - x[itr][0];
440  pz[itr][0] = zTmp - z[itr][0];
441  fGraphs[plane][itr]->GetPoint(fGraphs[plane][itr]->GetN()-1-fOffset,z[itr][1],x[itr][1]);
442  fGraphs[plane][itr]->GetPoint(fGraphs[plane][itr]->GetN()-2-fOffset,xTmp,zTmp);
443  px[itr][1] = xTmp - x[itr][1];
444  pz[itr][1] = zTmp - z[itr][1];
445  }
446  Int_t nLoops = 0;
447  while(1){
448  if(nLoops>10) break;
449  nLoops++;
450  Int_t nSplice = 0;
451 
452  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
453  if(spliced[itr1] != itr1) continue;
454  for(Int_t itr2 = itr1+1; itr2<fNtr[plane]; itr2++){
455  if(spliced[itr2] != itr2) continue;
456  Double_t dMin = 1e8;
457  Int_t iMin = -1, jMin = -1;
458  for(Int_t i = 0; i<2; i++){
459  for(Int_t j = 0; j<2; j++){
460  Double_t d = TMath::Abs(x[itr1][i]-x[itr2][j]) + TMath::Abs(z[itr1][i]-z[itr2][j]);
461  hDistTracks->Fill(d);
462  if(d<dMin){
463  dMin = d;
464  iMin = i;
465  jMin = j;
466  }
467  }
468  }
469  Double_t dx1 = x[itr1][1-iMin] - x[itr1][iMin];
470  Double_t dz1 = z[itr1][1-iMin] - z[itr1][iMin];
471  Double_t dx2 = x[itr2][1-jMin] - x[itr2][jMin];
472  Double_t dz2 = z[itr2][1-jMin] - z[itr2][jMin];
473 
474 
475  Double_t cross = dx1*dx2 + dz1*dz2;
476  Double_t d1 = dx1*dx1 + dz1*dz1;
477  Double_t d2 = dx2*dx2 + dz2*dz2;
478  Double_t sd1 = TMath::Sqrt(d1);
479  Double_t sd2 = TMath::Sqrt(d2);
480  if(d1*d2 == 0) continue;
481  Double_t costheta = cross/TMath::Sqrt(d1*d2);
482 
483 
484  if(gHarpoDebug>1)
485  Info("SpliceTracks","%d %d: (d=%g,theta=%g)",itr1,itr2,dMin,costheta);
486  // if(iMin==jMin) costheta = -costheta;
487  if(costheta > fThetaMin) continue;
488 
489 
490  hDistMinTracks->Fill(dMin);
491  hThetaTracks->Fill(1-costheta);
492  hThetaTracksNtr->Fill(dMin,1-costheta);
493  Double_t dForward1 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dx1 +
494  (z[itr1][iMin] - z[itr2][jMin])*dz1)/sd1;
495  Double_t dTransverse1 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dz1 -
496  (z[itr1][iMin] - z[itr2][jMin])*dx1)/sd1;
497  Double_t dForward2 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dx2 +
498  (z[itr1][iMin] - z[itr2][jMin])*dz2)/sd2;
499  Double_t dTransverse2 = TMath::Abs((x[itr1][iMin] - x[itr2][jMin])*dz2 -
500  (z[itr1][iMin] - z[itr2][jMin])*dx2)/sd2;
501 
502  if(gHarpoDebug>1)
503  Info("SpliceTracks","%d %d: (dF=%g,dT=%g,theta=%g)",itr1,itr2,dForward1,dTransverse1,costheta);
504 
505  if(fCanvas[plane]){
506  Double_t xRect[4] = {x[itr1][iMin] + ( fDminForward*dx1/sd1 - fDminTransverse*dz1/sd1),
507  x[itr1][iMin] + (- fDminForward*dx1/sd1 - fDminTransverse*dz1/sd1),
508  x[itr1][iMin] + (- fDminForward*dx1/sd1 + fDminTransverse*dz1/sd1),
509  x[itr1][iMin] + ( fDminForward*dx1/sd1 + fDminTransverse*dz1/sd1)};
510  Double_t zRect[4] = {z[itr1][iMin] + ( fDminForward*dz1/sd1 + fDminTransverse*dx1/sd1),
511  z[itr1][iMin] + (- fDminForward*dz1/sd1 + fDminTransverse*dx1/sd1),
512  z[itr1][iMin] + (- fDminForward*dz1/sd1 - fDminTransverse*dx1/sd1),
513  z[itr1][iMin] + ( fDminForward*dz1/sd1 - fDminTransverse*dx1/sd1)};
514 
515  std::cout << "----------------" << std::endl;
516  for(Int_t i = 0; i<4; i++) std::cout << xRect[i] << " ";
517  std::cout << std::endl;
518  for(Int_t i = 0; i<4; i++) std::cout << zRect[i] << " ";
519  std::cout << std::endl;
520  TPolyLine* rect = new TPolyLine(4,zRect,xRect,"FL");
521  rect->SetFillColor(kRed);
522  rect->SetLineWidth(2);
523  fCanvas[plane]->cd();
524  rect->Draw("L");
525  fHist[plane]->GetYaxis()->SetRangeUser(x[itr1][iMin]-20,x[itr1][iMin]+20);
526  fHist[plane]->GetXaxis()->SetRangeUser(z[itr1][iMin]-20,z[itr1][iMin]+20);
527  fCanvas[plane]->Update();
528  gSystem->Sleep(100);
529  }
530 
531 
532  if(nLoops<10)
533  hThetaDist[nLoops]->Fill(dMin,1-costheta);
534  if(((dForward1<fDminForward && dTransverse1<fDminTransverse) ||
535  (dForward2<fDminForward && dTransverse2<fDminTransverse)) &&
536  iMin >= 0){
537 
538 
539  if(fCanvas[plane]){
540  //a1 = 1; a2 = 1;
541  fCanvas[plane]->cd();
542  TArrow* arrow1 = new TArrow();
543  arrow1->SetLineColor(kBlue);
544  arrow1->DrawArrow(z[itr1][iMin],x[itr1][iMin],z[itr1][iMin] + dz1,x[itr1][iMin]+dx1);
545  TArrow* arrow2 = new TArrow();
546  arrow2->SetLineColor(kRed);
547  arrow2->DrawArrow(x[itr2][jMin],x[itr2][jMin],z[itr2][jMin]+dz2,x[itr2][jMin]+dx2);
548  fCanvas[plane]->Update();
549  gSystem->Sleep(1000);
550  fHist[plane]->GetXaxis()->SetRangeUser(90,420);
551  fHist[plane]->GetYaxis()->SetRangeUser(0,288);
552  }
553 
554 
555 
556  spliced[itr2] = spliced[itr1];
557  x[itr1][iMin] = x[itr2][1-jMin];
558  z[itr1][iMin] = z[itr2][1-jMin];
559  px[itr1][iMin] = px[itr2][1-jMin];
560  pz[itr1][iMin] = pz[itr2][1-jMin];
561  th[itr1][iMin] = th[itr2][1-jMin];
562  nSplice++;
563  if(gHarpoDebug>1)
564  Info("SpliceTrack","Splicing %d %d (dF=%g,dT=%g,theta=%g)",itr1,itr2,dForward1,dTransverse1,costheta);
565  continue;
566  }
567 
568  Double_t a1 = ((x[itr1][iMin] - x[itr2][jMin])*dz2 - (z[itr1][iMin] - z[itr2][jMin])*dx2)/(dz1*dx2 - dz2*dx1);
569  Double_t a2 = ((x[itr1][iMin] - x[itr2][jMin])*dz1 - (z[itr1][iMin] - z[itr2][jMin])*dx1)/(dz1*dx2 - dz2*dx1);
570 
571  if(gHarpoDebug>1)
572  Info("SpliceTracks","%d %d: %g/%g %g/%g",itr1,itr2,a1,sd1,a2,sd2);
573 
574 
575 
576  if( TMath::Abs(a1 - 0.5) > 1) continue;
577  if( TMath::Abs(a2 - 0.5) > 1) continue;
578  if( (( a1*sd1 < -40 ) || (a1*sd1 > sd1 + 40)) &&
579  (( a2*sd2 < -40 ) || (a2*sd2 > sd2 + 40)))
580  continue;
581 
582  if(fCanvas[plane]){
583  //a1 = 1; a2 = 1;
584  fCanvas[plane]->cd();
585  TArrow* arrow1 = new TArrow();
586  arrow1->SetLineColor(kBlue+3);
587  arrow1->DrawArrow(x[itr1][iMin],z[itr1][iMin],x[itr1][iMin] + a1*dx1,z[itr1][iMin]+a1*dz1);
588  TArrow* arrow2 = new TArrow();
589  arrow2->SetLineColor(kRed+3);
590  arrow2->DrawArrow(x[itr2][jMin],z[itr2][jMin],x[itr2][jMin]+a2*dx2,z[itr2][jMin]+a2*dz2);
591  fCanvas[plane]->Update();
592  gSystem->Sleep(1000);
593  fHist[plane]->GetXaxis()->SetRangeUser(90,420);
594  fHist[plane]->GetYaxis()->SetRangeUser(0,288);
595  }
596 
597  Double_t dMax = a1*sd1;
598  Int_t i0 = iMin, tr0 = itr1;
599  Int_t i1 = iMin, tr1 = itr1;
600  if((1-a1)*sd1 > dMax){
601  dMax = (1-a1)*sd1;
602  i1 = i0;
603  tr1 = tr0;
604  i0 = 1-iMin;
605  tr0 = itr1;
606  }
607  if(a2*sd2 > dMax){
608  dMax = a2*sd2;
609  i1 = i0;
610  tr1 = tr0;
611  i0 = jMin;
612  tr0 = itr2;
613  }
614  if((1-a2)*sd2 > dMax){
615  dMax = (1-a2)*sd2;
616  i1 = i0;
617  tr1 = tr0;
618  i0 = 1-jMin;
619  tr0 = itr2;
620  }
621 
622  spliced[itr2] = spliced[itr1];
623  x[itr1][0] = x[tr0][i0];
624  z[itr1][0] = z[tr0][i0];
625  px[itr1][0] = px[tr0][i0];
626  pz[itr1][0] = pz[tr0][i0];
627  th[itr1][0] = th[tr0][i0];
628  x[itr1][1] = x[tr1][i1];
629  z[itr1][1] = z[tr1][i1];
630  px[itr1][1] = px[tr1][i1];
631  pz[itr1][1] = pz[tr1][i1];
632  th[itr1][1] = th[tr1][i1];
633  nSplice++;
634  if(gHarpoDebug>1)
635  Info("SpliceTrack","** Splicing %d %d",itr1,itr2);
636  continue;
637  }
638  }
639  if(gHarpoDebug>1)
640  Info("SpliceTracks","** %d spliced (%d)",nSplice,nLoops);
641  if(nSplice == 0) break;
642  }
643 
644 
645 
646 
647 
648 
649  Int_t newTrIndex[fNtr[plane]];
650  Int_t nClTr[fNtr[plane]];
651  Double_t qTr[fNtr[plane]];
652  Int_t nSpliced = 0;
653  for(Int_t itr = 0; itr<fNtr[plane]; itr++){
654  nClTr[itr] = 0;
655  qTr[itr] = 0;
656  if(spliced[itr] != itr){
657  Int_t itrNew = spliced[itr];
658  while(spliced[itrNew] != itrNew)
659  itrNew = spliced[itrNew];
660  spliced[itr] = itrNew;
661  Int_t n = fGraphs[plane][itr]->GetN();
662  for(Int_t k = 0; k<n; k++){
663  Double_t xTmp, zTmp;
664  fGraphs[plane][itr]->GetPoint(k,xTmp,zTmp);
665  fGraphs[plane][itrNew]->SetPoint(fGraphs[plane][itrNew]->GetN(),xTmp,zTmp);
666  // fGraphs[plane][itr]->SetPoint(k,0,0);
667  }
668  fGraphs[plane][itr]->Set(0);
669  // fGraphs[plane][itr]->Delete();
670  }else{
671  nSpliced++;
672  }
673  newTrIndex[itr] = -1;//nSpliced;
674  }
675 
676 
677 
678  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
679  // Clusters in the event (both X and Y directions)
680  TClonesArray* clArray = reco->GetClustersArray();
681  Int_t nCl = clArray->GetEntries();
682  if(gHarpoDebug>1)
683  Info("Splice","ntot = %d",nCl);
684  // Int_t nClTot = 0, nClTot2 = 0;
685  Int_t nFinal = 0;
686  for(Int_t icl = 0; icl<nCl; icl++){
687  Int_t oldTrack = fTrack[icl];
688  if(oldTrack<0 || oldTrack>=fNtr[plane]) continue;
689  Int_t splicedTrack = spliced[oldTrack];
690  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
691  cluster->AddIdClusterTrack(splicedTrack);
692  nClTr[splicedTrack]++;
693  // nClNew[newTrackId]++;
694  }
695  GetQtracks(plane);
696  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
697  qTr[itr1] = fQcommon[itr1][itr1];
698  for(Int_t itr2 = itr1+1; itr2<fNtr[plane]; itr2++){
699  if(gHarpoDebug>1)
700  Info("SpliceTracks","%d => %d (%g %g %g)",itr2,spliced[itr1],fQcommon[itr1][itr1],fQcommon[itr2][itr2],fQcommon[itr1][itr2]);
701  if(!(fQcommon[itr1][itr2] > 0.4*fQcommon[itr2][itr2] || fQcommon[itr1][itr2] > 0.4*fQcommon[itr1][itr1] || fQcommon[itr1][itr2] > 20000))
702  continue;
703  if(gHarpoDebug>1)
704  Info("SpliceTracks","OK");
705  spliced[itr2] = spliced[itr1];
706  nClTr[spliced[itr1]] += nClTr[itr2];
707  qTr[spliced[itr1]] += fQcommon[itr2][itr2];
708  nSplice++;
709  }
710  }
711  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
712  if(itr1 == spliced[itr1] && nClTr[itr1]>fNclMin2 && qTr[itr1]>20000){
713  newTrIndex[itr1] = nFinal;
714  nFinal++;
715  }else
716  newTrIndex[itr1] = -1;
717 
718  }
719 
720  for(Int_t icl = 0; icl<nCl; icl++){
721  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
722  if(cluster->GetPlane() != plane) continue;
723  // Int_t nTr = cluster->GetNtr();
724  cluster->RemoveAllClusterTrack();
725  Int_t oldTrack = fTrack[icl];
726  if(oldTrack<0 || oldTrack>=fNtr[plane]) continue;
727  Int_t sp = spliced[oldTrack];
728  if(sp<0) continue;
729  if(newTrIndex[sp]<0) continue;
730  cluster->AddIdClusterTrack(newTrIndex[sp]);
731  }
732 
733  for(Int_t itr = 0; itr<fNtr[plane]; itr++){
734  if(gHarpoDebug>1)
735  Info("Splice","nCl(%d) = %d",itr,nClTr[itr]);
736 
737  if(spliced[itr] != itr) continue;
738  if(newTrIndex[itr] == -1) continue;
739  if(gHarpoDebug>1)
740  Info("Splice","nCl(%d) = %d / %d",itr,nClTr[itr],fGraphs[plane][itr]->GetN());
741  Double_t x0 = -1000, slope = -1000, chi2 = -1000;
742  if(fGraphs[plane][itr]->GetN()>10){
743  TLinearFitter* linfit = new TLinearFitter(1,"pol1","D");
744  linfit->AssignData(fGraphs[plane][itr]->GetN(), 1, fGraphs[plane][itr]->GetX(), fGraphs[plane][itr]->GetY());
745  linfit->EvalRobust();
746  x0 = linfit->GetParameter(0);
747  slope = linfit->GetParameter(1);
748  chi2 = linfit->GetChisquare();
749  linfit->Delete();
750  }
751 
752  Double_t qTrTmp = GetQtrack(itr,plane);
753  fQused += qTrTmp;
754 
755  if(gHarpoDebug>1)
756  Info("SpliceTracks","Qleft = %g / %g",fQtot-fQused,fQtot);
757 
758  //if(qTr<50000) continue;
759 
760  //Int_t ttz = 0, typeTrack = 0, IdTrackMatching = 0;
761  Int_t typeTrack = 0, IdTrackMatching = 0;
762  Double_t xSt = x[itr][0], zSt = z[itr][0], thSt = th[itr][0], xEnd = x[itr][1], zEnd = z[itr][1], thEnd = th[itr][1];
763  Double_t pxSt = px[itr][0], pzSt = pz[itr][0], pxEnd = px[itr][1], pzEnd = pz[itr][1];
764  if(gHarpoDebug>1)
765  Info("SpliceTracks","new Track: %d (%d cl, q=%g) => (%g, %g, %g) (%g, %g, %g)", newTrIndex[itr],nClTr[itr], qTrTmp, xSt, zSt, thSt, xEnd, zEnd, thEnd);
766  HarpoRecoKalmanTracks* tr = new HarpoRecoKalmanTracks(newTrIndex[itr],Int_t(qTrTmp),slope,x0,chi2,typeTrack,IdTrackMatching,nClTr[itr],plane,xSt,zSt,pxSt,pzSt,xEnd,zEnd,pxEnd,pzEnd);
767  // Info("SpliceTracks","Add graph %p",fGraphs[plane][itr]);
768  fEvt->GetRecoEvent()->AddTracks(tr);
769  fId[newTrIndex[itr]][plane] = itr;
770  tr->Delete();
771  }
772  if(gHarpoDebug>1)
773  Info("SpliceTracks","Done");
774 
775 }
776 
777 
778 
779 Double_t HarpoKalmanTracking::GetQtrack(Int_t itr, Int_t plane)
780 {
781  HarpoDccMap* map = fEvt->GetDccMap(plane);
782 
783  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
784  TClonesArray* clArray = reco->GetClustersArray();
785  Int_t nCl = clArray->GetEntries();
786  Double_t qTr = 0;
787  for(Int_t icl = 0; icl<nCl; icl++){
788  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
789  if(cluster->GetPlane() != plane) continue;
790  if(!cluster->CheckIdClusterTrack(itr)) continue;
791  Int_t type = cluster->GetType();
792  Int_t ind = cluster->GetIndex();
793  Int_t ind2 = cluster->GetIfirst();
794  Int_t width = cluster->GetWidth();
795  for(Int_t k = ind2; k< ind2 + width; k++){
796  Int_t i = -1, j = -1;
797  if(type == CCLUSTER){ i = ind; j = k;}
798  if(type == TCLUSTER){ i = k; j = ind;}
799  if(i<0 || i>=NCHAN) continue;
800  if(j<0 || j>=NADC) continue;
801  Double_t q = map->GetData(i,j);
802  // Double_t q = fMapTmp[i][j];
803  if(q<-500) continue;
804  ULong64_t test = 1 << itr;
805  //Info("GetQtrack","TEST = %llX | %llX => %llX, %llX",test,fMapTmp[i][j],fMapTmp[i][j]&test,fMapTmp[i][j]|test);
806  if(fMapTmp[i][j] & test) continue;
807  qTr += q;
808  // fMapTmp[i][j] = -1000;
809  fMapTmp[i][j] |= test;
810  }
811  }
812  if(gHarpoDebug>0)
813  Info("GetQtrack","Plane %d, Track %d: Q = %g/%g",plane,itr,qTr,fQtot);
814  return qTr;
815 }
816 
817 Double_t HarpoKalmanTracking::GetQtracks(Int_t plane)
818 {
819  HarpoDccMap* map = fEvt->GetDccMap(plane);
820 
821  HarpoRecoEvent* reco = fEvt->GetRecoEvent();
822  TClonesArray* clArray = reco->GetClustersArray();
823  Int_t nCl = clArray->GetEntries();
824  // Double_t qTr = 0;
825 
826  ULong_t used4Q[NCHAN][NADC];
827  for(Int_t i = 0; i<NCHAN; i++){
828  for(Int_t j = 0; j<NADC; j++){
829  // fMapTmp[i][j] = map->GetData(i,j);
830  used4Q[i][j] = 0;
831  }
832  }
833 
834  for(Int_t icl = 0; icl<nCl; icl++){
835  HarpoRecoClusters* cluster = (HarpoRecoClusters*)clArray->At(icl);
836  if(cluster->GetPlane() != plane) continue;
837  Int_t type = cluster->GetType();
838  Int_t ind = cluster->GetIndex();
839  Int_t ind2 = cluster->GetIfirst();
840  Int_t width = cluster->GetWidth();
841  Int_t nTr = cluster->GetNtr();
842  for(Int_t iTr = 0; iTr<nTr;iTr++){
843  Int_t itr = cluster->GetIdClusterTrack(iTr);
844  for(Int_t k = ind2; k< ind2 + width; k++){
845  Int_t i = -1, j = -1;
846  if(type == CCLUSTER){ i = ind; j = k;}
847  if(type == TCLUSTER){ i = k; j = ind;}
848  if(i<0 || i>=NCHAN) continue;
849  if(j<0 || j>=NADC) continue;
850  Double_t q = map->GetData(i,j);
851  // Double_t q = fMapTmp[i][j];
852  if(q<-500) continue;
853  ULong64_t test = 1 << itr;
854  //Info("GetQtrack","TEST = %llX | %llX => %llX, %llX",test,fMapTmp[i][j],fMapTmp[i][j]&test,fMapTmp[i][j]|test);
855  if(used4Q[i][j] & test) continue;
856  // qTr += q;
857  // fMapTmp[i][j] = -1000;
858  used4Q[i][j] |= test;
859  }
860  }
861  }
862 
863  // Double_t qTr[fNtr[plane]][fNtr[plane]];
864  for(Int_t i = 0; i<fNtr[plane]; i++){
865  for(Int_t j = 0; j<fNtr[plane]; j++){
866  fQcommon[i][j] = 0;
867  }
868  }
869  for(Int_t i = 0; i<NCHAN; i++){
870  for(Int_t j = 0; j<NADC; j++){
871  Double_t q = map->GetData(i,j);
872  if(q<-500) continue;
873 
874  for(Int_t itr1 = 0; itr1<fNtr[plane]; itr1++){
875  ULong64_t test1 = 1 << itr1;
876  if(!(used4Q[i][j] & test1)) continue;
877  for(Int_t itr2 = itr1; itr2<fNtr[plane]; itr2++){
878  ULong64_t test2 = 1 << itr2;
879  if(!(used4Q[i][j] & test2)) continue;
880  fQcommon[itr1][itr2] += q;
881  }
882  }
883  }
884  }
885 
886  if(gHarpoDebug>1){
887  for(Int_t i = 0; i<fNtr[plane]; i++){
888  for(Int_t j = 0; j<fNtr[plane]; j++){
889  std::cout << fQcommon[i][j] << " ";
890  }
891  std::cout << std::endl;
892  }
893  }
894 
895  // Info("GetQtrack","Plane %d, Track %d: Q = %g/%g",plane,itr,qTr,fQtot);
896  return 0;
897 }
898 
899 
901 {
903 
904  // Initialise histograms here
905 
906  fCanvas[0] = NULL;
907  fCanvas[1] = NULL;
908  fHist[0] = NULL;
909  fHist[1] = NULL;
910 
911  fNclMin = 10;
912  fNclMin2 = 15;
913  fMaxSlopeType = 1.5;
914  fResKalman = 0.2;
915  fScatKalman = 1;
916  fResFinder = 1;
917  fScatFinder = 0.5;
918 
919  fQmin = 200;
920  fWidthMinT = 2;
921  fWidthMinC = 5;
922  fWidthMax = 30;
923 
924 
925 
926 
927  // fNtr = 0;
928  hNcl = new TH1F("hNcl","hNcl",500,0,500);
929 
930  for(Int_t itr = 0; itr<NTRACK; itr++){
931  fGraphs[0][itr] = 0;
932  fGraphs[1][itr] = 0;
933  }
934 
935 
936  const Int_t nbinsDist = 200;
937  Double_t xminDist = 1;
938  Double_t xmaxDist = 10000;
939  Double_t logxminDist = TMath::Log(xminDist);
940  Double_t logxmaxDist = TMath::Log(xmaxDist);
941  Double_t binwidthDist = (logxmaxDist-logxminDist)/nbinsDist;
942  Double_t xbinsDist[nbinsDist+1];
943  xbinsDist[0] = xminDist;
944  for (Int_t i=1;i<=nbinsDist;i++)
945  xbinsDist[i] = TMath::Exp(logxminDist+i*binwidthDist);
946 
947 
948  const Int_t nbinsTheta = 100;
949  Double_t xminTheta = 1e-5;
950  Double_t xmaxTheta = 2;
951  Double_t logxminTheta = TMath::Log(xminTheta);
952  Double_t logxmaxTheta = TMath::Log(xmaxTheta);
953  Double_t binwidthTheta = (logxmaxTheta-logxminTheta)/nbinsTheta;
954  Double_t xbinsTheta[nbinsTheta+1];
955  xbinsTheta[0] = xminTheta;
956  for (Int_t i=1;i<=nbinsTheta;i++)
957  xbinsTheta[i] = TMath::Exp(logxminTheta+i*binwidthTheta);
958 
959  hDistTracks = new TH1F("hDistTracks","",nbinsDist,xbinsDist);
960  hDistMinTracks = new TH1F("hDistMinTracks","",nbinsDist,xbinsDist);
961  hThetaTracks = new TH1F("hThetaTracks","",nbinsTheta,xbinsTheta);//,100,-1,1);
962  // hThetaTracksNtr = new TH2F("hThetaTracksNtr","",NTRACK,0,NTRACK,100,-1,1);
963  hThetaTracksNtr = new TH2F("hThetaTracksNtr","",nbinsDist,xbinsDist,nbinsTheta,xbinsTheta);
964  for(Int_t i = 0; i<10; i++)
965  hThetaDist[i] = new TH2F(Form("hThetaDist%d",i),"",nbinsDist,xbinsDist,nbinsTheta,xbinsTheta);
966 
967 
968  const Int_t nbinsQ = 100;
969  Double_t xminQ = 1e-3;
970  Double_t xmaxQ = 1e3;
971  Double_t logxminQ = TMath::Log(xminQ);
972  Double_t logxmaxQ = TMath::Log(xmaxQ);
973  Double_t binwidthQ = (logxmaxQ-logxminQ)/nbinsQ;
974  Double_t xbinsQ[nbinsQ+1];
975  xbinsQ[0] = xminQ;
976  for (Int_t i=1;i<=nbinsQ;i++)
977  xbinsQ[i] = TMath::Exp(logxminQ+i*binwidthQ);
978  fHistQQtestVsDist = new TH2F("fHistQQtestVsDist","",nbinsQ,xbinsQ,nbinsDist,xbinsDist);
979  fHistQQmin = new TH2F("fHistQQmin","",nbinsQ,xbinsQ,nbinsDist,xbinsDist);
980 
981  }
982 
983 void HarpoKalmanTracking::Save(char * /* mode */)
984  {
985  OpenHistFile("recoKalmanTracking");
986  // Save histograms here
987  hDistTracks->Write();
988  hDistMinTracks->Write();
989  hThetaTracks->Write();
990  hThetaTracksNtr->Write();
991  for(Int_t i = 0; i<10; i++){
992  if(hThetaDist[i])
993  hThetaDist[i]->Write();
994  }
995  fHistQQtestVsDist->Write();
996  fHistQQmin->Write();
997 
998 
999  }
1000 
#define NTRACK
void print()
Ovreloaded medod whic do all job.
Int_t NCcl[NCHAN]
Definition: HarpoKalman.h:69
void SpliceTracks(Int_t plane)
HarpoRecoEvent * GetRecoEvent()
Definition: HarpoEvent.h:46
Double_t fScatKalman
Definition: HarpoKalman.h:128
Double_t fStartPointZ[NTRACK][10]
Definition: HarpoKalman.h:79
#define XPLANE
Definition: HarpoConfig.h:24
TH2F * hThetaDist[10]
Definition: HarpoKalman.h:107
#define NCHAN
Definition: HarpoDccMap.h:16
Int_t InitPlane(TClonesArray *clArray, Int_t plane)
Double_t GetMean()
Double_t fWidthMax
Definition: HarpoKalman.h:123
void AddIdClusterTrack(Int_t val)
Double_t fResKalman
Definition: HarpoKalman.h:127
Int_t fTrack[4000]
Definition: HarpoKalman.h:99
Int_t fId[NTRACK][2]
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
Definition: HarpoDccMap.cxx:84
TH1F * hThetaTracks
Definition: HarpoKalman.h:105
Double_t GetQtrack(Int_t itr, Int_t plane)
Int_t used[4000]
Definition: HarpoKalman.h:97
Double_t fQcommon[NTRACK][NTRACK]
TH2F * hThetaTracksNtr
Definition: HarpoKalman.h:106
Object containing the reconstruction information for one event (with array of HarpoRecoClusters Harpo...
Int_t AddTrack(TClonesArray *clArray, Int_t plane)
Double_t fMaxSlopeType
Definition: HarpoKalman.h:126
void AddTracks(HarpoRecoTracks *val)
Double_t fQmin
Definition: HarpoKalman.h:120
Int_t InitPlane(TClonesArray *clArray, Int_t plane)
Definition: HarpoKalman.cxx:54
TH1F * hDistMinTracks
Definition: HarpoKalman.h:104
TArrayI * FindNextClosest(TMatrixD X, TMatrixD C, Double_t angle, Int_t Ntr, TClonesArray *clArray, TArrayI *arr, Int_t ncl, Int_t plane, Int_t color, Int_t fill, Int_t smooth, Int_t q)
Int_t GetMapEdges(TClonesArray *clArray, Int_t plane, Int_t &iMin, Int_t &iMax, Int_t &jMin, Int_t &jMax)
Track finder with Kalman filter.
virtual void print()
FullEvent Header not scecific to the detectors The class is ....
TGraph * fGraphs[2][NTRACK]
TVirtualPad * fCanvas[2]
Definition: HarpoKalman.h:100
void RemoveAllClusterTrack()
Double_t GetQtracks(Int_t plane)
virtual void print()
Int_t fNclTrack[NTRACK]
Definition: HarpoKalman.h:74
HarpoRecoTracks object, Obtained with Kalman filter.
void Save(char *mode=NULL)
Cluster object, containing position, charge and quality information.
TFile * OpenHistFile(const char *ananame)
Double_t fScatFinder
Definition: HarpoKalman.h:130
static int type
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
ULong64_t fMapTmp[NCHAN][NADC]
Definition: HarpoKalman.h:93
Int_t NTcl[NADC]
Definition: HarpoKalman.h:71
Double_t fWidthMinC
Definition: HarpoKalman.h:122
Int_t GetIdClusterTrack()
#define NADC
Definition: HarpoDccMap.h:18
TH1F * hNcl
Redefine empty default.
#define KALMAN
Int_t fNtr[2]
Definition: HarpoKalman.h:72
Double_t fWidthMinT
Definition: HarpoKalman.h:121
#define TCLUSTER
Long64_t gHarpoDebug
Definition: HarpoDebug.cxx:9
Double_t fStartDirX[NTRACK][10]
Definition: HarpoKalman.h:81
Double_t fResFinder
Definition: HarpoKalman.h:129
TH1F * hDistTracks
Definition: HarpoKalman.h:103
Int_t fStartIndex[NTRACK]
Definition: HarpoKalman.h:87
HarpoEventHeader * GetHeader()
Definition: HarpoEvent.cxx:80
TArrayI * FindNext(TMatrixD X, TMatrixD C, Double_t angle, Int_t Ntr, TClonesArray *clArray, TArrayI *arr, Int_t ncl, Int_t plane, Int_t color, Int_t fill, Int_t smooth, Int_t q)
HarpoEvent * fEvt
Definition: HarpoAnalyse.h:70
ULong_t nEvents
Definition: HarpoAnalyse.h:75
Int_t Tcl[NADC][20]
Definition: HarpoKalman.h:70
Double_t fStartPointX[NTRACK][10]
Definition: HarpoKalman.h:78
const ULong_t gkNDetectors
Definition: HarpoDet.h:14
Double_t fStartDirZ[NTRACK][10]
Definition: HarpoKalman.h:82
#define CCLUSTER
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
Definition: HarpoEvent.cxx:108
void FindTrack(TClonesArray *clArray, Int_t icl0, Int_t icl1, TMatrixD Corig, Int_t plane, Int_t &color)
TClonesArray * GetClustersArray()
TH2F * fHist[2]
Definition: HarpoKalman.h:101
TH2F * fHistQQmin
Definition: HarpoKalman.h:110
void SetTrackType(Int_t val)
TH2F * fHistQQtestVsDist
Definition: HarpoKalman.h:109
Bool_t CheckIdClusterTrack(Int_t val)
Int_t Ccl[NCHAN][20]
Definition: HarpoKalman.h:68