--- calo/flight/CaloFranzini/src/Calib.cpp 2008/01/12 10:31:38 1.5 +++ calo/flight/CaloFranzini/src/Calib.cpp 2008/01/21 10:24:10 1.7 @@ -35,15 +35,15 @@ TMatrixD *matrix[17]; TMatrixD *nmat[17]; -TMatrixD *fqplane; -TMatrixD *fnqplane; -//TMatrixD *fqplane[17]; -//TMatrixD *fnqplane[17]; -//TMatrixF *fmatrix[17]; +//TMatrixD *fqplane; +//TMatrixD *fnqplane; +TMatrixD *fqplane[17]; +TMatrixD *fnqplane[17]; +TMatrixD *fmatrix[17]; //TMatrixF *fnmat[17]; //TMatrixD *fmatrix; //TMatrixD *fnmat; -TMatrixF *fmatrix; +//TMatrixF *fmatrix; //TMatrixF *fnmat; TMatrixF *fnmat[17]; //Int_t finmat[43][191]; @@ -251,7 +251,7 @@ }; // if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false; - if ( rigidity > 2.2 || rigidity < 1.5 ) return false; + // if ( rigidity > 2.2 || rigidity < 1.5 ) return false; // printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG); // return true; @@ -304,9 +304,10 @@ memset(rmean, 0, 17*sizeof(Float_t)); memset(ntot, 0, 17*sizeof(Int_t)); // memset(finmat, 0, 43*191*sizeof(Int_t)); +// Double_t tol = 1E-20; // -// for (Int_t i=0; i < 17 ; i++){ - for (Int_t i=3; i < 4 ; i++){ + for (Int_t i=0; i < 17 ; i++){ + // for (Int_t i=3; i < 4 ; i++){ if ( !FULL ){ matrix[i] = new TMatrixD(43,43); qplane[i] = new TArrayF(43); @@ -318,22 +319,31 @@ // fnmat = new TMatrixF(4128,4128); // fmatrix = new TMatrixF(8213,8213); // fnmat = new TMatrixF(8213,8213); - fmatrix = new TMatrixF(MDIM,MDIM); + // fmatrix = new TMatrixF(MDIM,MDIM); // fnmat = new TMatrixF(MDIM,MDIM); - fnmat[i] = new TMatrixF(43,191); -// cf->WriteFullMatrix(fmatrix, i); + // fmatrix[i] = new TMatrixF(1849,1849); + // fnmat[i] = new TMatrixF(43,43); + fmatrix[i] = new TMatrixD(1333,1333); +// fmatrix[i]->SetTol(tol); + fnmat[i] = new TMatrixF(43,31); + // cf->WriteFullMatrix(fmatrix, i); // cf->WriteFullNMatrix(fnmat, i); -// delete fmatrix; + // delete fmatrix; // delete fnmat; //fnmat[i] = new TMatrixI(8213,8213); } else { - 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) - fnqplane = new TMatrixD(43,191);// - // - cf->WriteFullMean(fqplane, i); - cf->WriteFullNMean(fnqplane, i); - delete fqplane; - delete fnqplane; + // 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) + // fnqplane = new TMatrixD(43,191);// + // fqplane[i] = new TMatrixD(43,43); // 43 planes x 43 "strip", where 43 = 50 + 14 + 5 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + [1] + ... + // fnqplane[i] = new TMatrixD(43,43);// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 + // + fqplane[i] = new TMatrixD(43,31); // 43 planes x 43 "strip", where 43 = 50 + 14 + 6 + 5 + 3 + 3 + 3 + 2 + 2 + 2 + 1 + 1 + 1 + 1 + 1 + [1] + 1 + 1 + 1 + 1 + ... + fnqplane[i] = new TMatrixD(43,31);// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ... + // + // cf->WriteFullMean(fqplane, i); + // cf->WriteFullNMean(fnqplane, i); + // delete fqplane; + // delete fnqplane; // }; }; @@ -394,8 +404,8 @@ // // FULL CALORIMETER // - fqplane = cf->LoadFullAverage(rbi); - fnqplane = cf->LoadFullNAverage(rbi); + // fqplane = cf->LoadFullAverage(rbi); + // fnqplane = cf->LoadFullNAverage(rbi); CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); // Int_t nplane = 0; @@ -409,7 +419,7 @@ Int_t mstrip = 0; // for (Int_t j=0; j<2; j++){ - for (Int_t i=0; i<21; i++){ + for (Int_t i=0; i<22; i++){ nplane = 1 - j + 2*i; if ( nplane > 37 ) nplane--; // @@ -417,10 +427,16 @@ // cd = 95 - cs; // + Int_t oldstr = -1; for (Int_t k=0; k<191; k++){ mstrip = cd + k; // if ( mstrip < (191-cs) ) (*fnqplane[rbi])[nplane][mstrip] += 1.; - if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.; + // if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.; + Int_t lstr = cf->ConvertStrip(mstrip); + if ( oldstr != lstr ){ + (*fnqplane[rbi])[nplane][lstr] += 1.; + oldstr = lstr; + }; }; }; }; @@ -439,17 +455,19 @@ // mstrip = cd + strip; // + Int_t lstr = cf->ConvertStrip(mstrip); // (*fqplane[rbi])[nplane][mstrip] += mip; - (*fqplane)[nplane][mstrip] += mip; + // (*fqplane)[nplane][mstrip] += mip; + (*fqplane[rbi])[nplane][lstr] += mip; // }; // - cf->WriteFullMean(fqplane, rbi); - cf->WriteFullNMean(fnqplane, rbi); - cf->UnLoadFullAverage(rbi); - cf->UnLoadFullNAverage(rbi); - delete fqplane; - delete fnqplane; + // cf->WriteFullMean(fqplane, rbi); + // cf->WriteFullNMean(fnqplane, rbi); + // cf->UnLoadFullAverage(rbi); + // cf->UnLoadFullNAverage(rbi); + // delete fqplane; + // delete fnqplane; // }; } @@ -476,40 +494,42 @@ } else { // for (Int_t i=0; iLoadFullAverage(i); - fnqplane = cf->LoadFullNAverage(i); + // fqplane = cf->LoadFullAverage(i); + // fnqplane = cf->LoadFullNAverage(i); if ( ntot[i] > 0 ) rmean[i] /= (Float_t)(ntot[i]); // for (Int_t j=0; j<43 ; j++){ - for (Int_t k=0; k<191; k++){ -// if ( (*fnqplane[i])[j][k] > 0 ){ -// (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k]; -// } else { -// (*fqplane[i])[j][k] = 0.; -// }; -// 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]); - if ( (*fnqplane)[j][k] > 0 ){ - if ( (*fqplane)[j][k] == 0. ) (*fqplane)[j][k] = 0.7; - (*fqplane)[j][k] /= (Float_t)(*fnqplane)[j][k]; + // for (Int_t k=0; k<191; k++){ + // for (Int_t k=0; k<43; k++){ + for (Int_t k=0; k<31; k++){ + // if ( (*fnqplane[i])[j][k] > 0 ){ + // (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k]; + // } else { + // (*fqplane[i])[j][k] = 0.; + // }; + // 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]); + if ( (*fnqplane[i])[j][k] > 0 ){ + if ( (*fqplane[i])[j][k] == 0. ) (*fqplane[i])[j][k] = 0.7; + (*fqplane[i])[j][k] /= (Float_t)(*fnqplane[i])[j][k]; } else { - (*fqplane)[j][k] = 0.; + (*fqplane[i])[j][k] = 0.; }; - 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]); + // 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]); }; }; - cf->WriteFullMean(fqplane, i); - cf->WriteFullNMean(fnqplane, i); - cf->UnLoadFullAverage(i); - cf->UnLoadFullNAverage(i); - delete fqplane; - delete fnqplane; + cf->WriteFullMean(fqplane[i], i); + cf->WriteFullNMean(fnqplane[i], i); + // cf->UnLoadFullAverage(i); + // cf->UnLoadFullNAverage(i); + // delete fqplane; + // delete fnqplane; }; // -// for (Int_t i=0; iWriteFullMean(fqplane[i], i); -// // -// }; + // for (Int_t i=0; iWriteFullMean(fqplane[i], i); + // // + // }; }; // cf->WriteNumBin(nbin); @@ -540,7 +560,6 @@ }; }; // -if ( rbi != 3 ) return; if ( erig < rig[0] ) return; if ( erig >= rig[nbin-1] ) return; // @@ -583,12 +602,13 @@ // // FULL CALORIMETER // - printf(" Retrieve matrix %i IEV %i \n",rbi,iev); + // if ( rbi != 3 ) return; + printf(" matrix %i IEV %i \n",rbi,iev); // fmatrix = cf->LoadFullMatrix(rbi); -// cf->LoadFullMatrix(rbi,fmatrix); + // cf->LoadFullMatrix(rbi,fmatrix); // fnmat = cf->LoadFullNMatrix(rbi); - printf(" done \n"); - printf(" start loop \n"); + // printf(" done \n"); + // printf(" start loop \n"); // CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); // @@ -600,16 +620,22 @@ // // Int_t mindgf = 48; // Int_t dgf = 143; + // Int_t mindgf = 0; //tutto + // Int_t dgf = 191; //tutto + // Int_t mindgf = 94; + // Int_t dgf = 96; + // Int_t mindgf = 84; + // Int_t dgf = 106; Int_t mindgf = 0; - Int_t dgf = 191; - // Int_t mindgf = 85; - // Int_t dgf = 115; + Int_t dgf = 43; Int_t cs = 0; Int_t cd = 0; Int_t mstrip = 0; // - Float_t mipv[43][191]; - memset(mipv,0,43*191*sizeof(Float_t)); + // Float_t mipv[43][43]; + // memset(mipv,0,43*43*sizeof(Float_t)); + Float_t mipv[43][31]; + memset(mipv,0,43*31*sizeof(Float_t)); // for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ // @@ -624,7 +650,8 @@ // mstrip = cd + strip; // - mipv[nplane][mstrip] = mip; + Int_t lstr = cf->ConvertStrip(mstrip); + mipv[nplane][lstr] += mip; // }; // @@ -659,10 +686,14 @@ // cd1 = 95 - cs1; // - mstrip1min = TMath::Max(mindgf,(cd1+0)); - mstrip1max = TMath::Min(dgf,(cd1+95)) + 1; + Int_t at1 = TMath::Max(0,(cd1+0)); + Int_t at2 = TMath::Min(190,(cd1+95)); + mstrip1min = cf->ConvertStrip(at1); + mstrip1max = cf->ConvertStrip(at2) + 1; + // mstrip1min = cf->ConvertStrip(TMath::Max(mindgf,(cd1+0))); + // mstrip1max = cf->ConvertStrip(TMath::Min(dgf,(cd1+95))) + 1; // - if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i \n",nplane1,mstrip1min,mstrip1max); + if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i mindgf %i dgf %i cd1 %i\n",nplane1,mstrip1min,mstrip1max,mindgf,dgf,cd1); // for (Int_t mstrip1=mstrip1min; mstrip1GetFullAverageAt(nplane1,mstrip1,erig,rbi); // - mi = (nplane1 * 191) + mstrip1; + // mi = (nplane1 * 191) + mstrip1; + // mi = (nplane1 * 43) + mstrip1; + mi = (nplane1 * 31) + mstrip1; // // if ( mstrip1 > mstrip1min ) break; // if ( mstrip1 > dgf ) break; @@ -695,18 +728,31 @@ // // mstrip2min = cd2 + 0; // mstrip2max = cd2 + 95; - mstrip2min = TMath::Max(mindgf,(cd2+0)); - mstrip2max = TMath::Min(dgf,(cd2+95)) + 1; + Int_t t1 = TMath::Max(0,(cd2+0)); + Int_t t2 = TMath::Min(190,(cd2+95)); + mstrip2min = cf->ConvertStrip(t1); + mstrip2max = cf->ConvertStrip(t2) + 1; // if ( nplane1 == 0 && nplane2 == 0 && mstrip1==mstrip1min ) printf(" mstrip2min %i mstrip2max %i \n",mstrip2min,mstrip2max); // for (Int_t mstrip2=mstrip2min; mstrip2GetFullAverageAt(nplane2,mstrip2,erig,rbi); // if ( mip2 != 0. ){ // - mj = (nplane2 * 191) + mstrip2; + // mj = (nplane2 * 191) + mstrip2; + // mj = (nplane2 * 43) + mstrip2; +// mj = (nplane2 * 31) + mstrip2; + Int_t sh = -15 + nplane1; + if ( sh > 15 ) sh -= 31*nplane1; + // + mj = (nplane2 * 31) + mstrip2 + sh; + // + if ( mj < 0 ) mj += 1333; + if ( mj >= 1333 ) mj -= 1333; +// printf(" mi %i mj %i sh %i \n",mi,mj,sh); + // // mj++; // // if ( mstrip2 > mstrip2min ) break; @@ -715,7 +761,8 @@ // if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){ // (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig)); // (*fnmat[rbi])[mi][mj] += 1.; - (*fmatrix)[mi][mj] += (mip1 * mip2); // giusto + (*fmatrix[rbi])[mi][mj] += (mip1 * mip2); // giusto + // (*fmatrix)[mi][mj] += (mip1 * mip2) * 1000000.; toto++; // (*fmatrix)[mi][mj] += 1.; // cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi); @@ -731,18 +778,18 @@ // printf(" toto = %i \n",toto); printf("\n done \n"); - printf(" write matrix \n"); -// cf->WriteFullMatrix(fmatrix, rbi); + // printf(" write matrix \n"); + // cf->WriteFullMatrix(fmatrix, rbi); // cf->WriteFullNMatrix(fnmat, rbi); - printf(" done \n"); - printf(" unload matrix \n"); - // cf->UnLoadFullMatrix(rbi); + // printf(" done \n"); + // printf(" unload matrix \n"); + // cf->UnLoadFullMatrix(rbi); // cf->UnLoadFullNMatrix(rbi); - printf(" done \n"); - printf(" delete matrix \n"); -// delete fmatrix; + // printf(" done \n"); + // printf(" delete matrix \n"); + // delete fmatrix; // delete fnmat; - printf(" done \n"); + // printf(" done \n"); }; } @@ -798,112 +845,253 @@ // // FULL // - // for (Int_t i=0; iLoadFullMatrix(i); + // fmatrix = cf->LoadFullMatrix(i); // fnmat = cf->LoadFullNMatrix(i); // -// for (Int_t ii=0; ii 0. ){ -// // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; -// // } else { -// // (*fmatrix[i])[ii][j] = 0.; -// // }; -// if ( (*fnmat)[ii][j] > 0. ){ -// (*fmatrix)[ii][j] /= (*fnmat)[ii][j]; -// } else { -// (*fmatrix)[ii][j] = 0.; -// }; -// }; -// }; -// + // for (Int_t ii=0; ii 0. ){ + // // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; + // // } else { + // // (*fmatrix[i])[ii][j] = 0.; + // // }; + // if ( (*fnmat)[ii][j] > 0. ){ + // (*fmatrix)[ii][j] /= (*fnmat)[ii][j]; + // } else { + // (*fmatrix)[ii][j] = 0.; + // }; + // }; + // }; + // +// TMatrixD *mymat3 = new TMatrixD(129,129); +// TMatrixD *mymat5 = new TMatrixD(215,215); +// TMatrixD *mymat7 = new TMatrixD(301,301); +// TMatrixD *mymat9 = new TMatrixD(387,387); +// TMatrixD *mymat11 = new TMatrixD(473,473); +// TMatrixD *mymat17 = new TMatrixD(731,731); + // TMatrixF *mymat = new TMatrixF(129,129); + // TMatrixF *mymat = new TMatrixF(989,989); Int_t i1 = -1; Int_t j1 = -1; + // int mi,mj; Int_t nonzero = 0; Int_t nonzero1 = 0; for (Int_t ii=0; ii<43; ii++){ - for (Int_t j=0; j<191; j++){ -// if ( (*fnmat[i])[ii][j] > 0. ){ -// (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; -// } else { -// (*fmatrix[i])[ii][j] = 0.; -// }; - i1 = (ii * 191) + j; + // for (Int_t j=0; j<191; j++){ + // for (Int_t j=0; j<43; j++){ + for (Int_t j=0; j<31; j++){ + // if ( (*fnmat[i])[ii][j] > 0. ){ + // (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; + // } else { + // (*fmatrix[i])[ii][j] = 0.; + // }; + // i1 = (ii * 191) + j; + // i1 = (ii * 43) + j; + i1 = (ii * 31) + j; // j1 = -1; for (Int_t iij=0; iij<43; iij++){ - for (Int_t jj=0; jj<191; jj++){ + // for (Int_t jj=0; jj<191; jj++){ + // for (Int_t jj=0; jj<43; jj++){ + for (Int_t jj=0; jj<31; jj++){ + // + // j1 = (iij * 191) + jj; + // j1 = (iij * 43) + jj; + Int_t sh = -15 + ii; + if ( sh > 15 ) sh -= 31*ii; + // + j1 = (iij * 31) + jj + sh; // - j1 = (iij * 191) + jj; + if ( j1 < 0 ) j1 += 1333; + if ( j1 >= 1333 ) j1 -= 1333; + +// j1 = (iij * 31) + jj; // j1++; // if ( finmat[ii][j] > 0 ){ // (*fmatrix)[i1][j1] /= finmat[ii][j]; - if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1] == (*fmatrix)[i1][j1]) ){ - (*fmatrix)[i1][j1] = 0.; + if ( (*fnmat[i])[ii][j] == 0. || (*fmatrix[i])[i1][j1] == 0. || !((*fmatrix[i])[i1][j1] == (*fmatrix[i])[i1][j1]) ){ + (*fmatrix[i])[i1][j1] = 1.; } else { - (*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j]; + (*fmatrix[i])[i1][j1] /= (*fnmat[i])[ii][j]; nonzero++; if ( i1 == 0 ) nonzero1++; }; + // +// if ( j>=7 && j <=23 && jj >=7 && jj<=23 ){ +// Int_t mi17 = (ii*3) + j -7; +// Int_t mj17 = (iij*3) + jj -7; +// (*mymat17)[mi17][mj17] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=10 && j <=20 && jj >=10 && jj<=20 ){ +// Int_t mi11 = (ii*3) + j -10; +// Int_t mj11 = (iij*3) + jj -10; +// (*mymat11)[mi11][mj11] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=11 && j <=19 && jj >=11 && jj<=19 ){ +// Int_t mi9 = (ii*3) + j -11; +// Int_t mj9 = (iij*3) + jj -11; +// (*mymat9)[mi9][mj9] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=12 && j <=18 && jj >=12 && jj<=18 ){ +// Int_t mi7 = (ii*3) + j -12; +// Int_t mj7 = (iij*3) + jj -12; +// (*mymat7)[mi7][mj7] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=13 && j <=17 && jj >=13 && jj<=17 ){ +// Int_t mi5 = (ii*3) + j -13; +// Int_t mj5 = (iij*3) + jj -13; +// (*mymat5)[mi5][mj5] = (*fmatrix[i])[i1][j1]; +// }; +// if ( j>=14 && j <=16 && jj >=14 && jj<=16 ){ +// Int_t mi3 = (ii*3) + j -14; +// Int_t mj3 = (iij*3) + jj -14; +// (*mymat3)[mi3][mj3] = (*fmatrix[i])[i1][j1]; +// }; + + + // if ( j>=94 && j <=96 && jj >=94 && jj<=96 ){ + // mi = (ii*3) + j -94; + // mj = (iij*3) + jj -94; + // (*mymat)[mi][mj] = (*fmatrix)[i1][j1]; + // }; + + + // if ( j>=84 && j <=106 && jj >=84 && jj<=106 ){ + // mi = (ii*3) + j -84; + // mj = (iij*3) + jj -84; + // (*mymat)[mi][mj] = (*fmatrix)[i1][j1]; + // }; + }; }; }; }; // printf(" Matrix has %i non-zero elements \n",nonzero); - printf(" Matrix has %i non-zero elements on the first row\n",nonzero1); + // printf(" Matrix has %i non-zero elements on the first row\n",nonzero1); // -// Bool_t BAD = false; -// for (Int_t ii=0; ii<43; ii++){ -// for (Int_t j=0; j<191; j++){ -// // -// i1 = (ii * 191) + j; -// // -// for (Int_t iij=0; iij<43; iij++){ -// for (Int_t jj=0; jj<191; jj++){ -// // -// j1 = (iij * 191) + jj; -// // -// // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); -// if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){ -// printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); -// printf(" che schifo! \n"); -// BAD = true; -// }; -// // -// }; -// }; -// }; -// }; -// // -// if ( BAD ) printf(" questa matrice fa cagare \n"); + // Bool_t BAD = false; + // for (Int_t ii=0; ii<43; ii++){ + // for (Int_t j=0; j<191; j++){ + // // + // i1 = (ii * 191) + j; + // // + // for (Int_t iij=0; iij<43; iij++){ + // for (Int_t jj=0; jj<191; jj++){ + // // + // j1 = (iij * 191) + jj; + // // + // // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); + // if ( (*fmatrix)[i1][j1] == 0. || !((*fmatrix)[i1][j1]==(*fmatrix)[i1][j1]) ){ + // printf(" ROW %i COLUMN %i VALUE %f \n",i1,j1,(*fmatrix)[i1][j1]); + // printf(" che schifo! \n"); + // BAD = true; + // }; + // // + // }; + // }; + // }; + // }; + // // + // if ( BAD ) printf(" questa matrice fa cagare \n"); // // - // cf->WriteFullMatrix(fmatrix[i],i); - cf->WriteFullMatrix(fmatrix, i); + cf->WriteFullMatrix(fmatrix[i],i); + // cf->WriteFullMatrix(fmatrix, i); // cf->WriteFullNMatrix(fnmat, i); cf->WriteFullNMatrix(fnmat[i], i); // - // if ( fmatrix[i]->Determinant() == 0. ){ - if ( fmatrix->Determinant() == 0. ){ + // TDecompSVD svd(*fmatrix[i]); + // Bool_t ok = svd.Decompose(); + // + Double_t zero = (Double_t)0.0; + // + if ( fmatrix[i]->Determinant() == zero ){ + //if ( fmatrix->Determinant() == 0. ){ printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); } else { - Double_t det = 0.; - // TMatrixF invmatrix = (TMatrixF)(fmatrix[i]->Invert(&det)); - // printf(" Bin %i determinant is %f \n",i,det); - // cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); - TMatrixF invmatrix = (TMatrixF)(fmatrix->Invert(&det)); - printf(" Bin %i determinant is %f \n",i,det); - cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); + // }; + // if ( i == 3 ){ + // if ( ok ){ + // Double_t tol = 1E-20; + // TDecompSVD svd((*fmatrix)[i],tol); + // svd.Decompose(); + // TMatrixD svdInv = svd.Invert(); + // svdInv.Print("svdInv"); + // cout << "condition: " << svd.Condition() << endl; + // cf->WriteInvertedFullMatrix((TMatrixD)svdInv,999); + + Double_t det = 0.; + TMatrixD invmatrix = (TMatrixD)(fmatrix[i]->Invert(&det)); + printf(" Bin %i determinant is %f \n",i,det); + cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,i); }; + +// if ( mymat3->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat3->Invert(&det)); +// printf(" Mymat3 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1103); +// }; +// cf->WriteFullMatrix(mymat3, 103); +// if ( mymat5->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat5->Invert(&det)); +// printf(" Mymat5 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1105); +// }; +// cf->WriteFullMatrix(mymat5, 105); +// if ( mymat7->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat7->Invert(&det)); +// printf(" Mymat7 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1107); +// }; +// cf->WriteFullMatrix(mymat7, 107); +// if ( mymat9->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat9->Invert(&det)); +// printf(" Mymat3 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1109); +// }; +// cf->WriteFullMatrix(mymat9, 109); +// if ( mymat11->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat11->Invert(&det)); +// printf(" Mymat11 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1111); +// }; +// cf->WriteFullMatrix(mymat11, 111); +// if ( mymat17->Determinant() == 0. ){ +// printf(" ERROR: the matrix at bin %i is singular, determinant = 0., it cannot be inverted! \n",i); +// } else { +// Double_t det = 0.; +// TMatrixD invmatrix = (TMatrixD)(mymat17->Invert(&det)); +// printf(" Mymat3 determinant is %f \n",det); +// cf->WriteInvertedFullMatrix((TMatrixD)invmatrix,1117); +// }; +// cf->WriteFullMatrix(mymat17, 117); + + // - cf->UnLoadFullMatrix(i); + // cf->UnLoadFullMatrix(i); // cf->UnLoadFullNMatrix(i); - delete fmatrix; + // delete fmatrix; // delete fnmat; // };