25 #include "TRootEmbeddedCanvas.h"
26 #include "TGListBox.h"
55 if (plane != NULL )plane->
print();
68 Int_t nV = vArray->GetEntries();
78 for(Int_t j = 0; j<3; j++){
84 for(Int_t iV = 0; iV<nV; iV++){
92 Double_t mean[7][nV*nV];
93 Double_t sigma[7][nV*nV];
94 Double_t cov[7][nV*nV];
96 Int_t iX[nV*nV], iY[nV*nV];
97 for(Int_t i = 0; i<nV*nV; i++){
100 for(Int_t j = 0; j<7; j++){
106 for(Int_t iVx = 0; iVx<nV; iVx++){
110 for(Int_t iVy = 0; iVy<nV; iVy++){
114 if(TMath::Abs(zVx - zVy)>20)
continue;
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]);
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]);
152 Int_t* index =
new Int_t[nPair];
154 TMath::Sort(nPair,cov[0],index,kTRUE);
167 Double_t norm[2] = {0,0};
168 for(Int_t i = 0; i<nMCtr; 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]);
177 for(Int_t k = 0; k<3; k++){
178 psim[0][k] /= norm[0];
179 psim[1][k] /= norm[1];
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;
214 Info(
"process",
"Keeping pair %d - %d",iX[index[i]],iY[index[i]]);
218 Bool_t test = cov[5][index[i]] < cov[6][index[i]];
223 chi2 = cov[3][index[i]]+cov[4][index[i]];
225 Info(
"process",
"1->2 2->1");
228 Info(
"process",
"1->1 2->2");
230 chi2 = cov[1][index[i]]+cov[2][index[i]];
233 if(cov[5][index[i]] == cov[6][index[i]])
234 j = Int_t(gRandom->Rndm()+0.5);
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;
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;
271 if( (psim[0][0]-psim[1][0])*(psim[0][1]-psim[1][1])*(pxX[0]-pxX[1])*(pxY[0]-pxY[1])>0 )
277 j = Int_t(gRandom->Rndm()+0.5);
280 if(pzX[0]>0 && pzY[j]>0){
281 p1[0] = pxX[0]/pzX[0];
282 p1[1] = pxY[j]/pzY[j];
284 }
else if(pzX[0]<0 && pzY[j]<0){
285 p1[0] = -pxX[0]/pzX[0];
286 p1[1] = -pxY[j]/pzY[j];
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];
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];
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];
312 norm1 = TMath::Sqrt(norm1);
313 norm2 = TMath::Sqrt(norm2);
314 for(Int_t k = 0; k<3; k++){
319 data[0] = 0.5*(zVx+zVy);
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];
340 if( (psim[0][0]-psim[1][0])*(psim[0][1]-psim[1][1])*(pxX[0]-pxX[1])*(pxY[0]-pxY[1])>0 )
344 if((psim[0][0]-psim[1][0])*(psim[0][1]-psim[1][1])*(pxX[0]-pxX[1])*(pxY[0]-pxY[1])==0)
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;
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);
366 Int_t nbins = hX->GetXaxis()->GetNbins();
367 if(hY->GetXaxis()->GetNbins() != nbins)
return -1;
370 TArrayD* arr =
new TArrayD(nbins);
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);
388 if(qXold>0) dqX = (qX - qXold)/(qX+qXold);
389 if(qYold>0) dqY = (qY - qYold)/(qY+qYold);
395 if(dqY==-10000)
continue;
396 if(dqX==-10000)
continue;
398 if(hoverlapX->GetBinContent(bin)>0) w *= 0.1;
399 if(hoverlapY->GetBinContent(bin)>0) w *= 0.1;
402 arr->AddAt(qX/qY,bin);
405 sigma += w*qX/qY*qX/qY;
437 cov /= TMath::Sqrt(sigX*sigY);
440 sigma = sigma/W - mean*mean;
442 Info(
"CompareHist",
"n = %d, W = %g, mean = %g, sigma = %g",n, W, mean, sigma);
463 Int_t nbins = hX1->GetXaxis()->GetNbins();
464 if(hY1->GetXaxis()->GetNbins() != nbins)
return -1;
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];
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]);
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;
498 if(hoverlapX->GetBinContent(bin)>0) w *= 0.1;
499 if(hoverlapY->GetBinContent(bin)>0) w *= 0.1;
505 cov += w*dqX[i]*dqY[i];
507 sigX += w*dqX[i]*dqX[i];
509 sigY += w*dqY[i]*dqY[i];
524 cov /= TMath::Sqrt(sigX*sigY);
549 for(Int_t i=0;i<
NADC;i++){
550 for(Int_t j=0;j<
NCHAN;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){
564 hoverlap->Fill(i,500);
567 h2->Fill(i,q*(1.-r));
586 Int_t nbins = h1->GetXaxis()->GetNbins();
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;
604 if(dq1==-10000)
continue;
605 if(dq2==-10000)
continue;
622 Long64_t perfectMatching = 0;
623 if ( !
gHConfig->
Lookup(
"matching.matchingMethod",perfectMatching) )
624 Info(
"Init",
"Use normal matching");
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;
639 for(Int_t j = 0; j<3; j++){
640 hQvT[i][j] =
new TH1F(Form(
"hQvT_%d_%d",i,j),
"",512,0,512);
642 hOverlap[i] =
new TH1F(Form(
"hOverlap_%d",i),
"",512,0,512);
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);
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");
672 TCanvas* c = ecTab->GetCanvas();
673 c->GetPad(1)->Delete();
674 c->GetPad(2)->Delete();
679 Int_t nV = vArray->GetEntries();
680 Int_t nVp[2] = {0,0};
682 for(Int_t i = 0; i<nV; i++){
685 iV[plane][nVp[plane]] = i;
704 TLatex* latex =
new TLatex();
705 latex->SetTextFont(132);
706 latex->SetTextAlign(12);
707 latex->SetTextSize(0.08);
712 Int_t nV3d = v3dArray->GetEntries();
713 c->GetPad(1)->Divide(4,nV3d+1);
714 c->GetPad(2)->Divide(4,nV3d+1);
723 infobox->NewEntry(Form(
"%d vertexes, %d 3D vertexes",nV,nV3d));
724 for(Int_t k = 0; k<nV3d; k++){
726 infobox->NewEntry(Form(
"%d: (%g,%g,%g) 1:(%g,%g,%g) 2:(%g,%g,%g)",
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]);
776 Bool_t save = kFALSE;
777 if(save) outfile =
new TFile(Form(
"exampleMatching_%d.root",k),
"recreate");
781 if(save) hCorrel11->Write(
"hCorrel11");
787 if(save) hCorrel22->Write(
"hCorrel22");
788 hCorrel22->Add(hCorrel11);
791 if(save) hCorrel12->Write(
"hCorrel12");
797 if(save) hCorrel21->Write(
"hCorrel21");
798 hCorrel21->Add(hCorrel12);
801 latex->DrawLatex(500,5000,Form(
"%g",TMath::Max(cov11,cov22)));
803 latex->DrawLatex(500,5000,Form(
"%g",TMath::Max(cov21,cov12)));
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++){
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);
830 for(Int_t i=0;i<
NADC;i++){
831 for(Int_t j=0;j<
NCHAN;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));
852 h1->SetLineColor(kGreen);
853 h2->SetLineColor(kMagenta);
855 h2->DrawCopy(
"box same");
856 latex->DrawLatex(50,200,Form(
"%d",iV[plane][k]));
882 if(save) h1->Write();
883 if(save) h2->Write();
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");
896 if(save) outfile->Close();
898 hQvT[iVx][1]->SetLineColor(kRed);
899 hQvT[iVx][2]->SetLineColor(kRed);
900 hQvT[iVy][1]->SetLineColor(kBlue);
901 hQvT[iVy][2]->SetLineColor(kBlue);
909 hQvT[iVy][1]->DrawCopy(
"same");
914 hQvT[iVy][2]->DrawCopy(
"same");
919 hQvT[iVy][2]->DrawCopy(
"same");
924 hQvT[iVy][1]->DrawCopy(
"same");
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();
1016 main->SetCleanup(kDeepCleanup);
1018 TGVerticalFrame* frame =
new TGVerticalFrame(main,xsize,ysize,kVerticalFrame);
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);
1032 TGLabel* fAnalysisLabel =
new TGLabel(frame,
"Template analysis");
1035 TGHorizontalFrame* nBinsFrame =
new TGHorizontalFrame(frame,xsize,ysize2,kHorizontalFrame);
1036 TGLabel* nBinsLabel =
new TGLabel(nBinsFrame,
"fNbins");
1040 main->AddFrame(frame,fLayout4);
1041 frame->AddFrame(fAnalysisLabel,fLayout3);
1042 frame->AddFrame(nBinsFrame,fLayout3);
1043 nBinsFrame->AddFrame(nBinsLabel,fLayout1);
1050 TGTextButton* setConf =
new TGTextButton(frame,
"Save Config",
id);
1051 setConf->Associate(fMain);
1053 frame->AddFrame(setConf,fLayout3);
1055 main->MapSubwindows();
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()
Double_t Covariance(TH1F *hX1, TH1F *hY1, TH1F *hX2, TH1F *hY2, TH1F *hoverlapX, TH1F *hoverlapY)
void DisplayAnalysis(TRootEmbeddedCanvas *ecTab, TGListBox *infobox)
void Save(char *mode=NULL)
Double_t GetPz1Vertex3D()
Double_t GetPy2Vertex3D()
void print()
Overloaded method which do all job.
Double_t GetData(Int_t i, Int_t j)
Set/Get Data Cell.
int main(int argc, char **argv)
Double_t GetThetaVertex()
Bool_t isExist(ULong_t det)
Detecror date exist //! Number of Real Detectors.
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)
TClonesArray * GetVertexArray()
A class store HARPO row DCC event data and header. End provide access metods to the row data...
Double_t GetPx1Vertex3D()
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...
HarpoDetEvent * GetDetEvent(Long_t plane=XDCC)
TH1F * hOverlap[kMaxNvertex]
Double_t GetPy1Vertex3D()
Data from Keller temperuture and pressure sensors.
TH1F * hQvT[kMaxNvertex][3]
void ResetVertex3DArray()
TVirtualPad * MakeNiceHisto(TH1 *hist, TVirtualPad *c1, const char *opt, Bool_t copy)
TpcSimMCTrack * GetTrack(Int_t i)
TpcSimMCEvent * GetMCEvent()
Bool_t Lookup(const char *path, Bool_t &val)
Lookup function for scalar values.
TClonesArray * GetVertex3DArray()
static const Int_t kMaxNvertex
HarpoEventHeader * GetHeader()
Double_t GetPx2Vertex3D()
Double_t GetPz2Vertex3D()
const ULong_t gkNDetectors
HarpoDccMap * GetDccMap(Long_t plane=XDCC)
R__EXTERN HarpoDetSet * gHDetSet
3D vertex object, containing position, angle and associated 2D vertexes, and quality info ...