| 25 |
Float_t rig[18]; |
Float_t rig[18]; |
| 26 |
Float_t rmean[17]; |
Float_t rmean[17]; |
| 27 |
Int_t ntot[17]; |
Int_t ntot[17]; |
| 28 |
|
Int_t MDIM = 8213; |
| 29 |
|
//Int_t MDIM = 4128; |
| 30 |
//Float_t qqplane[17][43]; |
//Float_t qqplane[17][43]; |
| 31 |
//Int_t nnqplane[17][43]; |
//Int_t nnqplane[17][43]; |
| 32 |
|
|
| 35 |
TMatrixD *matrix[17]; |
TMatrixD *matrix[17]; |
| 36 |
TMatrixD *nmat[17]; |
TMatrixD *nmat[17]; |
| 37 |
|
|
| 38 |
TMatrixD *fqplane[17]; |
TMatrixD *fqplane; |
| 39 |
TMatrixD *fnqplane[17]; |
TMatrixD *fnqplane; |
| 40 |
TMatrixF *fmatrix[17]; |
//TMatrixD *fqplane[17]; |
| 41 |
|
//TMatrixD *fnqplane[17]; |
| 42 |
|
//TMatrixF *fmatrix[17]; |
| 43 |
|
//TMatrixF *fnmat[17]; |
| 44 |
|
//TMatrixD *fmatrix; |
| 45 |
|
//TMatrixD *fnmat; |
| 46 |
|
TMatrixF *fmatrix; |
| 47 |
|
//TMatrixF *fnmat; |
| 48 |
TMatrixF *fnmat[17]; |
TMatrixF *fnmat[17]; |
| 49 |
|
//Int_t finmat[43][191]; |
| 50 |
|
|
| 51 |
//=============================================================================== |
//=============================================================================== |
| 52 |
bool Select( PamLevel2* event ){ |
bool Select( PamLevel2* event ){ |
| 250 |
if ( view == 1 && (plane >0 || strip > ti0) ) break; |
if ( view == 1 && (plane >0 || strip > ti0) ) break; |
| 251 |
}; |
}; |
| 252 |
// if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false; |
// if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false; |
| 253 |
|
|
| 254 |
|
// if ( rigidity > 2.2 || rigidity < 1.5 ) return false; |
| 255 |
|
// printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG); |
| 256 |
// |
// |
| 257 |
return true; |
return true; |
| 258 |
} |
} |
| 303 |
// |
// |
| 304 |
memset(rmean, 0, 17*sizeof(Float_t)); |
memset(rmean, 0, 17*sizeof(Float_t)); |
| 305 |
memset(ntot, 0, 17*sizeof(Int_t)); |
memset(ntot, 0, 17*sizeof(Int_t)); |
| 306 |
|
// memset(finmat, 0, 43*191*sizeof(Int_t)); |
| 307 |
// |
// |
| 308 |
for (Int_t i=0; i < 17 ; i++){ |
for (Int_t i=0; i < 17 ; i++){ |
| 309 |
|
//for (Int_t i=3; i < 4 ; i++){ |
| 310 |
if ( !FULL ){ |
if ( !FULL ){ |
| 311 |
matrix[i] = new TMatrixD(43,43); |
matrix[i] = new TMatrixD(43,43); |
| 312 |
qplane[i] = new TArrayF(43); |
qplane[i] = new TArrayF(43); |
| 314 |
nmat[i] = new TMatrixD(43,43); |
nmat[i] = new TMatrixD(43,43); |
| 315 |
} else { |
} else { |
| 316 |
if ( MATRIX ){ |
if ( MATRIX ){ |
| 317 |
fmatrix[i] = new TMatrixF(4128,4128); |
// fmatrix = new TMatrixF(4128,4128); |
| 318 |
fnmat[i] = new TMatrixF(4128,4128); |
// fnmat = new TMatrixF(4128,4128); |
| 319 |
|
// fmatrix = new TMatrixF(8213,8213); |
| 320 |
|
// fnmat = new TMatrixF(8213,8213); |
| 321 |
|
fmatrix = new TMatrixF(MDIM,MDIM); |
| 322 |
|
// fnmat = new TMatrixF(MDIM,MDIM); |
| 323 |
|
fnmat[i] = new TMatrixF(43,191); |
| 324 |
|
// cf->WriteFullMatrix(fmatrix, i); |
| 325 |
|
// cf->WriteFullNMatrix(fnmat, i); |
| 326 |
|
// delete fmatrix; |
| 327 |
|
// delete fnmat; |
| 328 |
//fnmat[i] = new TMatrixI(8213,8213); |
//fnmat[i] = new TMatrixI(8213,8213); |
| 329 |
} else { |
} else { |
| 330 |
fqplane[i] = new TMatrixD(43,191); // 43 planes x 191 strip (= 1 + 95 x 2, one strip is the one transversed by the track that could be on the extreme right or left) |
fqplane = new TMatrixD(43,191); // 43 planes x 191 strip (= 1 + 95 x 2, one strip is the one transversed by the track that could be on the extreme right or left) |
| 331 |
fnqplane[i] = new TMatrixD(43,191);// |
fnqplane = new TMatrixD(43,191);// |
| 332 |
|
// |
| 333 |
|
cf->WriteFullMean(fqplane, i); |
| 334 |
|
cf->WriteFullNMean(fnqplane, i); |
| 335 |
|
delete fqplane; |
| 336 |
|
delete fnqplane; |
| 337 |
|
// |
| 338 |
}; |
}; |
| 339 |
}; |
}; |
| 340 |
}; |
}; |
| 394 |
// |
// |
| 395 |
// FULL CALORIMETER |
// FULL CALORIMETER |
| 396 |
// |
// |
| 397 |
|
fqplane = cf->LoadFullAverage(rbi); |
| 398 |
|
fnqplane = cf->LoadFullNAverage(rbi); |
| 399 |
CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); |
CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); |
| 400 |
// |
// |
| 401 |
Int_t nplane = 0; |
Int_t nplane = 0; |
| 419 |
// |
// |
| 420 |
for (Int_t k=0; k<191; k++){ |
for (Int_t k=0; k<191; k++){ |
| 421 |
mstrip = cd + k; |
mstrip = cd + k; |
| 422 |
if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.; |
// if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.; |
| 423 |
|
if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.; |
| 424 |
}; |
}; |
| 425 |
}; |
}; |
| 426 |
}; |
}; |
| 439 |
// |
// |
| 440 |
mstrip = cd + strip; |
mstrip = cd + strip; |
| 441 |
// |
// |
| 442 |
(*fqplane[rbi])[nplane][mstrip] += mip; |
// (*fqplane[rbi])[nplane][mstrip] += mip; |
| 443 |
|
(*fqplane)[nplane][mstrip] += mip; |
| 444 |
// |
// |
| 445 |
}; |
}; |
| 446 |
|
// |
| 447 |
|
cf->WriteFullMean(fqplane, rbi); |
| 448 |
|
cf->WriteFullNMean(fnqplane, rbi); |
| 449 |
|
cf->UnLoadFullAverage(rbi); |
| 450 |
|
cf->UnLoadFullNAverage(rbi); |
| 451 |
|
delete fqplane; |
| 452 |
|
delete fnqplane; |
| 453 |
// |
// |
| 454 |
}; |
}; |
| 455 |
} |
} |
| 474 |
// |
// |
| 475 |
}; |
}; |
| 476 |
} else { |
} else { |
| 477 |
|
// |
| 478 |
for (Int_t i=0; i<nbin-1; i++){ |
for (Int_t i=0; i<nbin-1; i++){ |
| 479 |
|
fqplane = cf->LoadFullAverage(i); |
| 480 |
|
fnqplane = cf->LoadFullNAverage(i); |
| 481 |
if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]); |
if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]); |
| 482 |
// |
// |
| 483 |
for (Int_t j=0; j<43 ; j++){ |
for (Int_t j=0; j<43 ; j++){ |
| 484 |
for (Int_t k=0; k<191; k++){ |
for (Int_t k=0; k<191; k++){ |
| 485 |
if ( (*fnqplane[i])[j][k] > 0 ){ |
// if ( (*fnqplane[i])[j][k] > 0 ){ |
| 486 |
(*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k]; |
// (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k]; |
| 487 |
|
// } else { |
| 488 |
|
// (*fqplane[i])[j][k] = 0.; |
| 489 |
|
// }; |
| 490 |
|
// printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane[i])[j][k],(*fnqplane[i])[j][k]); |
| 491 |
|
if ( (*fnqplane)[j][k] > 0 ){ |
| 492 |
|
if ( (*fqplane)[j][k] == 0. ) (*fqplane)[j][k] = 0.7; |
| 493 |
|
(*fqplane)[j][k] /= (Float_t)(*fnqplane)[j][k]; |
| 494 |
} else { |
} else { |
| 495 |
(*fqplane[i])[j][k] = 0.; |
(*fqplane)[j][k] = 0.; |
| 496 |
}; |
}; |
| 497 |
printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane[i])[j][k],(*fnqplane[i])[j][k]); |
printf(" BIN %i plane %i strip %i average energy %f qplane %f nqplane %f \n",i,j,k,rmean[i],(*fqplane)[j][k],(*fnqplane)[j][k]); |
| 498 |
}; |
}; |
| 499 |
}; |
}; |
| 500 |
|
cf->WriteFullMean(fqplane, i); |
| 501 |
|
cf->WriteFullNMean(fnqplane, i); |
| 502 |
|
cf->UnLoadFullAverage(i); |
| 503 |
|
cf->UnLoadFullNAverage(i); |
| 504 |
|
delete fqplane; |
| 505 |
|
delete fnqplane; |
| 506 |
}; |
}; |
| 507 |
for (Int_t i=0; i<nbin-1; i++){ |
// |
| 508 |
// |
// for (Int_t i=0; i<nbin-1; i++){ |
| 509 |
cf->WriteFullMean(fqplane[i], i); |
// // |
| 510 |
// |
// cf->WriteFullMean(fqplane[i], i); |
| 511 |
}; |
// // |
| 512 |
|
// }; |
| 513 |
}; |
}; |
| 514 |
// |
// |
| 515 |
cf->WriteNumBin(nbin); |
cf->WriteNumBin(nbin); |
| 540 |
}; |
}; |
| 541 |
}; |
}; |
| 542 |
// |
// |
| 543 |
|
if ( rbi != 3 ) return; |
| 544 |
if ( erig < rig[0] ) return; |
if ( erig < rig[0] ) return; |
| 545 |
if ( erig >= rig[nbin-1] ) return; |
if ( erig >= rig[nbin-1] ) return; |
| 546 |
// |
// |
| 583 |
// |
// |
| 584 |
// FULL CALORIMETER |
// FULL CALORIMETER |
| 585 |
// |
// |
| 586 |
|
printf(" Retrieve matrix %i IEV %i \n",rbi,iev); |
| 587 |
|
// fmatrix = cf->LoadFullMatrix(rbi); |
| 588 |
|
// cf->LoadFullMatrix(rbi,fmatrix); |
| 589 |
|
// fnmat = cf->LoadFullNMatrix(rbi); |
| 590 |
|
printf(" done \n"); |
| 591 |
|
printf(" start loop \n"); |
| 592 |
|
// |
| 593 |
CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); |
CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); |
| 594 |
// |
// |
| 595 |
Int_t nplane = 0; |
Int_t nplane = 0; |
| 598 |
Int_t strip = 0; |
Int_t strip = 0; |
| 599 |
Float_t mip = 0.; |
Float_t mip = 0.; |
| 600 |
// |
// |
| 601 |
Int_t mindgf = 48; |
// Int_t mindgf = 48; |
| 602 |
Int_t dgf = 143; |
// Int_t dgf = 143; |
| 603 |
|
Int_t mindgf = 0; |
| 604 |
|
Int_t dgf = 191; |
| 605 |
|
// Int_t mindgf = 85; |
| 606 |
|
// Int_t dgf = 115; |
| 607 |
Int_t cs = 0; |
Int_t cs = 0; |
| 608 |
Int_t cd = 0; |
Int_t cd = 0; |
| 609 |
Int_t mstrip = 0; |
Int_t mstrip = 0; |
| 628 |
// |
// |
| 629 |
}; |
}; |
| 630 |
// |
// |
| 631 |
Float_t mip1; |
Float_t mip1 = 1.; |
| 632 |
Int_t cs1; |
Int_t cs1; |
| 633 |
Int_t cd1; |
Int_t cd1; |
| 634 |
Float_t mip2; |
Float_t mip2 = 1.; |
| 635 |
Int_t cs2; |
Int_t cs2; |
| 636 |
Int_t cd2; |
Int_t cd2; |
| 637 |
Int_t mi = -1; |
Int_t mi = -1; |
| 647 |
Int_t mstrip2min = 0; |
Int_t mstrip2min = 0; |
| 648 |
Int_t mstrip2max = 0; |
Int_t mstrip2max = 0; |
| 649 |
// |
// |
| 650 |
|
Int_t toto = 0; |
| 651 |
|
// |
| 652 |
for (Int_t nplane1=0; nplane1<43; nplane1++){ |
for (Int_t nplane1=0; nplane1<43; nplane1++){ |
| 653 |
for (Int_t mstrip1=0; mstrip1<191; mstrip1++){ |
if ( nplane1 >= 37 ) nn1 = nplane1 + 1; |
| 654 |
|
vi1 = 1; |
| 655 |
|
if ( nn1%2 ) vi1 = 0; |
| 656 |
|
pl1 = (nn1 - 1 + vi1)/2; |
| 657 |
|
// |
| 658 |
|
cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1 |
| 659 |
|
// |
| 660 |
|
cd1 = 95 - cs1; |
| 661 |
|
// |
| 662 |
|
mstrip1min = TMath::Max(mindgf,(cd1+0)); |
| 663 |
|
mstrip1max = TMath::Min(dgf,(cd1+95)) + 1; |
| 664 |
|
// |
| 665 |
|
if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i \n",nplane1,mstrip1min,mstrip1max); |
| 666 |
|
// |
| 667 |
|
for (Int_t mstrip1=mstrip1min; mstrip1<mstrip1max; mstrip1++){ |
| 668 |
|
// printf(".\n"); |
| 669 |
|
// |
| 670 |
mj = -1; |
mj = -1; |
| 671 |
// |
// |
| 672 |
if ( nplane1 >= 37 ) nn1 = nplane1 + 1; |
mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi); |
| 673 |
vi1 = 1; |
// |
| 674 |
if ( nn1%2 ) vi1 = 0; |
mi = (nplane1 * 191) + mstrip1; |
| 675 |
pl1 = (nn1 - 1 + vi1)/2; |
// |
| 676 |
// |
// if ( mstrip1 > mstrip1min ) break; |
| 677 |
cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1 |
// if ( mstrip1 > dgf ) break; |
| 678 |
// |
// if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){ |
| 679 |
cd1 = 95 - cs1; |
// |
| 680 |
// |
// finmat[nplane1][mstrip1]++; |
| 681 |
mstrip1min = cd1 + 0; |
(*fnmat[rbi])[nplane1][mstrip1] += 1.; |
| 682 |
mstrip1max = cd1 + 95; |
// |
| 683 |
mip1 = mipv[nplane1][mstrip1]; |
if ( mip1 != 0. ){ |
| 684 |
// |
// |
|
if ( mstrip1 >= mindgf && mstrip1 <= dgf ){ |
|
|
mi++; |
|
|
// |
|
| 685 |
for (Int_t nplane2=0; nplane2<43; nplane2++){ |
for (Int_t nplane2=0; nplane2<43; nplane2++){ |
| 686 |
for (Int_t mstrip2=0; mstrip2<191; mstrip2++){ |
// |
| 687 |
// |
if ( nplane2 >= 37 ) nn2 = nplane2 + 1; |
| 688 |
if ( nplane2 >= 37 ) nn2 = nplane2 + 1; |
vi2 = 1; |
| 689 |
vi2 = 1; |
if ( nn2%2 ) vi2 = 0; |
| 690 |
if ( nn2%2 ) vi2 = 0; |
pl1 = (nn2 - 1 + vi2)/2; |
| 691 |
pl1 = (nn2 - 1 + vi2)/2; |
// |
| 692 |
|
cs2 = ct->tibar[pl2][vi2] - 1; |
| 693 |
|
// |
| 694 |
|
cd2 = 95 - cs2; |
| 695 |
|
// |
| 696 |
|
// mstrip2min = cd2 + 0; |
| 697 |
|
// mstrip2max = cd2 + 95; |
| 698 |
|
mstrip2min = TMath::Max(mindgf,(cd2+0)); |
| 699 |
|
mstrip2max = TMath::Min(dgf,(cd2+95)) + 1; |
| 700 |
|
// |
| 701 |
|
if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max); |
| 702 |
|
// |
| 703 |
|
for (Int_t mstrip2=mstrip2min; mstrip2<mstrip2max; mstrip2++){ |
| 704 |
// |
// |
| 705 |
cs2 = ct->tibar[pl2][vi2] - 1; |
mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi); |
| 706 |
// |
// |
| 707 |
cd2 = 95 - cs2; |
if ( mip2 != 0. ){ |
| 708 |
// |
// |
| 709 |
mstrip2min = cd2 + 0; |
mj = (nplane2 * 191) + mstrip2; |
| 710 |
mstrip2max = cd2 + 95; |
// mj++; |
| 711 |
// |
// |
| 712 |
mip2 = mipv[nplane2][mstrip2]; |
// if ( mstrip2 > mstrip2min ) break; |
| 713 |
// |
// if ( mstrip2 > dgf ) break; |
| 714 |
if ( mstrip2 >= mindgf && mstrip2 <= dgf ){ |
// if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){ |
| 715 |
mj++; |
// if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){ |
| 716 |
(*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig)); |
// (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig)); |
| 717 |
if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){ |
// (*fnmat[rbi])[mi][mj] += 1.; |
| 718 |
(*fnmat[rbi])[mi][mj] += 1.; |
(*fmatrix)[mi][mj] += (mip1 * mip2); // giusto |
| 719 |
}; |
toto++; |
| 720 |
|
// (*fmatrix)[mi][mj] += 1.; |
| 721 |
|
// cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi); |
| 722 |
|
// cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi); |
| 723 |
|
// (*fnmat)[mi][mj] += 1.; |
| 724 |
|
// }; |
| 725 |
}; |
}; |
| 726 |
}; |
}; |
| 727 |
}; |
}; |
| 729 |
}; |
}; |
| 730 |
}; |
}; |
| 731 |
// |
// |
| 732 |
|
printf(" toto = %i \n",toto); |
| 733 |
|
printf("\n done \n"); |
| 734 |
|
printf(" write matrix \n"); |
| 735 |
|
// cf->WriteFullMatrix(fmatrix, rbi); |
| 736 |
|
// cf->WriteFullNMatrix(fnmat, rbi); |
| 737 |
|
printf(" done \n"); |
| 738 |
|
printf(" unload matrix \n"); |
| 739 |
|
// cf->UnLoadFullMatrix(rbi); |
| 740 |
|
// cf->UnLoadFullNMatrix(rbi); |
| 741 |
|
printf(" done \n"); |
| 742 |
|
printf(" delete matrix \n"); |
| 743 |
|
// delete fmatrix; |
| 744 |
|
// delete fnmat; |
| 745 |
|
printf(" done \n"); |
| 746 |
}; |
}; |
| 747 |
} |
} |
| 748 |
|
|
| 798 |
// |
// |
| 799 |
// FULL |
// FULL |
| 800 |
// |
// |
| 801 |
for (Int_t i=0; i<nbin-1; i++){ |
// for (Int_t i=0; i<nbin-1; i++){ |
| 802 |
|
// for (Int_t i=3; i<5; i++){ |
| 803 |
|
for (Int_t i=3; i<4; i++){ |
| 804 |
// |
// |
| 805 |
// determine the average matrix |
// determine the average matrix |
| 806 |
// |
// |
| 807 |
for (Int_t ii=0; ii<4171; ii++){ |
// fmatrix = cf->LoadFullMatrix(i); |
| 808 |
for (Int_t j=0; j<4171; j++){ |
// fnmat = cf->LoadFullNMatrix(i); |
| 809 |
if ( (*fnmat[i])[ii][j] > 0. ){ |
// |
| 810 |
(*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; |
// for (Int_t ii=0; ii<MDIM; ii++){ |
| 811 |
} else { |
// for (Int_t j=0; j<MDIM; j++){ |
| 812 |
(*fmatrix[i])[ii][j] = 0.; |
// // if ( (*fnmat[i])[ii][j] > 0. ){ |
| 813 |
|
// // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; |
| 814 |
|
// // } else { |
| 815 |
|
// // (*fmatrix[i])[ii][j] = 0.; |
| 816 |
|
// // }; |
| 817 |
|
// if ( (*fnmat)[ii][j] > 0. ){ |
| 818 |
|
// (*fmatrix)[ii][j] /= (*fnmat)[ii][j]; |
| 819 |
|
// } else { |
| 820 |
|
// (*fmatrix)[ii][j] = 0.; |
| 821 |
|
// }; |
| 822 |
|
// }; |
| 823 |
|
// }; |
| 824 |
|
// |
| 825 |
|
Int_t i1 = -1; |
| 826 |
|
Int_t j1 = -1; |
| 827 |
|
Int_t nonzero = 0; |
| 828 |
|
Int_t nonzero1 = 0; |
| 829 |
|
for (Int_t ii=0; ii<43; ii++){ |
| 830 |
|
for (Int_t j=0; j<191; j++){ |
| 831 |
|
// if ( (*fnmat[i])[ii][j] > 0. ){ |
| 832 |
|
// (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; |
| 833 |
|
// } else { |
| 834 |
|
// (*fmatrix[i])[ii][j] = 0.; |
| 835 |
|
// }; |
| 836 |
|
i1 = (ii * 191) + j; |
| 837 |
|
// j1 = -1; |
| 838 |
|
for (Int_t iij=0; iij<43; iij++){ |
| 839 |
|
for (Int_t jj=0; jj<191; jj++){ |
| 840 |
|
// |
| 841 |
|
j1 = (iij * 191) + jj; |
| 842 |
|
// j1++; |
| 843 |
|
// if ( finmat[ii][j] > 0 ){ |
| 844 |
|
// (*fmatrix)[i1][j1] /= finmat[ii][j]; |
| 845 |
|
if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1] == (*fmatrix)[i1][j1]) ){ |
| 846 |
|
(*fmatrix)[i1][j1] = 0.; |
| 847 |
|
} else { |
| 848 |
|
(*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j]; |
| 849 |
|
nonzero++; |
| 850 |
|
if ( i1 == 0 ) nonzero1++; |
| 851 |
|
}; |
| 852 |
|
}; |
| 853 |
}; |
}; |
| 854 |
}; |
}; |
| 855 |
}; |
}; |
| 856 |
// |
// |
| 857 |
cf->WriteFullMatrix(fmatrix[i],i); |
printf(" Matrix has %i non-zero elements \n",nonzero); |
| 858 |
|
printf(" Matrix has %i non-zero elements on the first row\n",nonzero1); |
| 859 |
|
// |
| 860 |
|
// Bool_t BAD = false; |
| 861 |
|
// for (Int_t ii=0; ii<43; ii++){ |
| 862 |
|
// for (Int_t j=0; j<191; j++){ |
| 863 |
|
// // |
| 864 |
|
// i1 = (ii * 191) + j; |
| 865 |
|
// // |
| 866 |
|
// for (Int_t iij=0; iij<43; iij++){ |
| 867 |
|
// for (Int_t jj=0; jj<191; jj++){ |
| 868 |
|
// // |
| 869 |
|
// j1 = (iij * 191) + jj; |
| 870 |
|
// // |
| 871 |
|
// // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); |
| 872 |
|
// if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){ |
| 873 |
|
// printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); |
| 874 |
|
// printf(" che schifo! \n"); |
| 875 |
|
// BAD = true; |
| 876 |
|
// }; |
| 877 |
|
// // |
| 878 |
|
// }; |
| 879 |
|
// }; |
| 880 |
|
// }; |
| 881 |
|
// }; |
| 882 |
|
// // |
| 883 |
|
// if ( BAD ) printf(" questa matrice fa cagare \n"); |
| 884 |
|
// |
| 885 |
|
// |
| 886 |
|
// cf->WriteFullMatrix(fmatrix[i],i); |
| 887 |
|
cf->WriteFullMatrix(fmatrix, i); |
| 888 |
|
// cf->WriteFullNMatrix(fnmat, i); |
| 889 |
|
cf->WriteFullNMatrix(fnmat[i], i); |
| 890 |
// |
// |
| 891 |
if ( fmatrix[i]->Determinant() == 0. ){ |
// if ( fmatrix[i]->Determinant() == 0. ){ |
| 892 |
|
if ( fmatrix->Determinant() == 0. ){ |
| 893 |
printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); |
printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); |
| 894 |
} else { |
} else { |
| 895 |
Double_t det = 0.; |
Double_t det = 0.; |
| 896 |
TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det)); |
// TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det)); |
| 897 |
|
// printf(" Bin %i determinant is %f \n",i,det); |
| 898 |
|
// cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); |
| 899 |
|
TMatrixF invmatrix = (TMatrixF)(fmatrix->Invert(&det)); |
| 900 |
printf(" Bin %i determinant is %f \n",i,det); |
printf(" Bin %i determinant is %f \n",i,det); |
| 901 |
cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); |
cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); |
| 902 |
}; |
}; |
| 903 |
|
// |
| 904 |
|
cf->UnLoadFullMatrix(i); |
| 905 |
|
// cf->UnLoadFullNMatrix(i); |
| 906 |
|
delete fmatrix; |
| 907 |
|
// delete fnmat; |
| 908 |
|
// |
| 909 |
}; |
}; |
| 910 |
}; |
}; |
| 911 |
// |
// |