--- calo/flight/CaloFranzini/src/Calib.cpp 2008/01/11 15:27:13 1.4 +++ calo/flight/CaloFranzini/src/Calib.cpp 2008/01/25 15:09:07 1.8 @@ -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]; @@ -58,6 +58,7 @@ PamTrack *track = event->GetTrack(0); if(!track)return false; + if ( SIMU ) return true; //------------------------------------------------------------------ // tracker pre-selection //------------------------------------------------------------------ @@ -264,16 +265,22 @@ // // //=============================================================== -void CreateHistos( PamLevel2* event , TString file){ +void CreateHistos( PamLevel2* event , TString file){ cf = new CaloFranzini(event); // + if ( FULL ){ + cf->CalculateLongTZeta(false); + cf->CalculateFullTZeta(); + }; + // if ( MATRIX ){ cf->UpdateMatrixFile(file.Data()); - cf->LoadBin(); if ( !FULL ){ + cf->LoadBin(); cf->LoadLong(); } else { + cf->LoadBin(true); cf->LoadFull(); }; } else { @@ -304,9 +311,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=3; i < 4 ; i++){ if ( !FULL ){ matrix[i] = new TMatrixD(43,43); qplane[i] = new TArrayF(43); @@ -318,22 +326,35 @@ // 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); + // fnmat[i] = new TMatrixF(43,31); + fmatrix[i] = new TMatrixD(473,473); + fnmat[i] = new TMatrixF(43,11); +// fmatrix[i]->SetTol(tol); + // 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 ... + fqplane[i] = new TMatrixD(43,11); // 43 planes x 11 "strip", where 43 = 84 + 6 + 3 + 1 + 1 +[1]+ 1 + 1 + 3 + 6 + 84 + fnqplane[i] = new TMatrixD(43,11);// 0 1 2 3 4 5 6 7 8 9 10 + // + // cf->WriteFullMean(fqplane, i); + // cf->WriteFullNMean(fnqplane, i); + // delete fqplane; + // delete fnqplane; // }; }; @@ -394,8 +415,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 +430,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 +438,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 +466,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 +505,43 @@ } 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++){ + for (Int_t k=0; k<11; 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); @@ -518,9 +550,15 @@ cf->WriteRigBin(rigbin); // TArrayF *rmeanbin = new TArrayF(17, rmean); - TFile *file = cf->GetFile(); - file->cd(); - file->WriteObject(&(*rmeanbin), "binrigmean"); + if ( FULL ){ + TFile *file = cf->GetFullFile(); + file->cd(); + file->WriteObject(&(*rmeanbin), "binrigmean"); + } else { + TFile *file = cf->GetFile(); + file->cd(); + file->WriteObject(&(*rmeanbin), "binrigmean"); + }; // // } @@ -540,7 +578,6 @@ }; }; // -if ( rbi != 3 ) return; if ( erig < rig[0] ) return; if ( erig >= rig[nbin-1] ) return; // @@ -583,12 +620,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 +638,25 @@ // // 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 dgf = 11; 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)); + Float_t mipv[43][11]; + memset(mipv,0,43*11*sizeof(Float_t)); // for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ // @@ -624,7 +671,8 @@ // mstrip = cd + strip; // - mipv[nplane][mstrip] = mip; + Int_t lstr = cf->ConvertStrip(mstrip); + mipv[nplane][lstr] += mip; // }; // @@ -649,6 +697,17 @@ // Int_t toto = 0; // + // + // + Int_t therigb = rbi; + // + if ( erig < rig[0] ){ + therigb = 0; + }; + if ( erig >= rig[nbin-2] ){ + therigb = nbin - 3; + }; + // for (Int_t nplane1=0; nplane1<43; nplane1++){ if ( nplane1 >= 37 ) nn1 = nplane1 + 1; vi1 = 1; @@ -659,19 +718,26 @@ // 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); + mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,therigb); // - mi = (nplane1 * 191) + mstrip1; + // mi = (nplane1 * 191) + mstrip1; + // mi = (nplane1 * 43) + mstrip1; + // mi = (nplane1 * 31) + mstrip1; + mi = (nplane1 * 11) + mstrip1; // // if ( mstrip1 > mstrip1min ) break; // if ( mstrip1 > dgf ) break; @@ -695,18 +761,46 @@ // // 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); + // + mip2 = mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig,therigb); // 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; +// +// +// +// +// Int_t sh = -5 + nplane1; +// if ( sh > 5 ) sh -= 11*nplane1; +// // +// mj = (nplane2 * 11) + mstrip2 + sh; +// // +// if ( mj < 0 ) mj += 473; +// if ( mj >= 473 ) mj -= 473; +// // + // + mj = (nplane2 * 11) + mstrip2; + // +// printf(" mi %i mj %i sh %i \n",mi,mj,sh); + // // mj++; // // if ( mstrip2 > mstrip2min ) break; @@ -715,7 +809,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 +826,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 +893,268 @@ // // 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++){ + for (Int_t j=0; j<11; 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; + i1 = (ii * 11) + 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++){ + for (Int_t jj=0; jj<11; jj++){ // - j1 = (iij * 191) + 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; +// // +// if ( j1 < 0 ) j1 += 1333; +// if ( j1 >= 1333 ) j1 -= 1333; + // + // + // +// Int_t sh = -5 + ii; +// if ( sh > 5 ) sh -= 11*ii; +// // +// j1 = (iij * 11) + jj + sh; +// // +// if ( j1 < 0 ) j1 += 473; +// if ( j1 >= 473 ) j1 -= 473; +// + j1 = (iij * 11) + 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); // cf->WriteFullNMatrix(fnmat, i); - cf->WriteFullNMatrix(fnmat[i], i); + //cf->WriteFullMatrix(fmatrix[i],i);// OK + // cf->WriteFullNMatrix(fnmat[i], i); //OK + // + // TDecompSVD svd(*fmatrix[i]); + // Bool_t ok = svd.Decompose(); // - // if ( fmatrix[i]->Determinant() == 0. ){ - if ( fmatrix->Determinant() == 0. ){ + 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); + cf->WriteFullMatrix(fmatrix[i],i);// OK + cf->WriteFullNMatrix(fnmat[i], i); //OK } 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; // };