--- calo/flight/CaloFranzini/src/Calib.cpp 2008/01/03 10:02:26 1.3 +++ calo/flight/CaloFranzini/src/Calib.cpp 2008/01/15 12:41:38 1.6 @@ -25,6 +25,8 @@ Float_t rig[18]; Float_t rmean[17]; Int_t ntot[17]; +Int_t MDIM = 8213; +//Int_t MDIM = 4128; //Float_t qqplane[17][43]; //Int_t nnqplane[17][43]; @@ -33,11 +35,18 @@ TMatrixD *matrix[17]; TMatrixD *nmat[17]; -TMatrixD *fqplane[17]; -TMatrixD *fnqplane[17]; -TMatrixF *fmatrix[17]; +TMatrixD *fqplane; +TMatrixD *fnqplane; +//TMatrixD *fqplane[17]; +//TMatrixD *fnqplane[17]; +//TMatrixF *fmatrix[17]; +//TMatrixF *fnmat[17]; +//TMatrixD *fmatrix; +//TMatrixD *fnmat; +TMatrixF *fmatrix; +//TMatrixF *fnmat; TMatrixF *fnmat[17]; - +//Int_t finmat[43][191]; //=============================================================================== bool Select( PamLevel2* event ){ @@ -241,6 +250,9 @@ if ( view == 1 && (plane >0 || strip > ti0) ) break; }; // if ( event->GetCaloLevel1()->qtotpl(0) > 7. ) return false; + + if ( rigidity > 2.2 || rigidity < 1.5 ) return false; + // printf(" rig %f CRIG %i SRIG %i \n",rigidity,CRIG,SRIG); // return true; } @@ -291,8 +303,10 @@ // memset(rmean, 0, 17*sizeof(Float_t)); memset(ntot, 0, 17*sizeof(Int_t)); + // memset(finmat, 0, 43*191*sizeof(Int_t)); // - for (Int_t i=0; i < 17 ; 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); @@ -300,12 +314,27 @@ nmat[i] = new TMatrixD(43,43); } else { if ( MATRIX ){ - fmatrix[i] = new TMatrixF(4128,4128); - fnmat[i] = new TMatrixF(4128,4128); + // fmatrix = new TMatrixF(4128,4128); + // fnmat = new TMatrixF(4128,4128); + // fmatrix = new TMatrixF(8213,8213); + // fnmat = new TMatrixF(8213,8213); + fmatrix = new TMatrixF(MDIM,MDIM); + // fnmat = new TMatrixF(MDIM,MDIM); + fnmat[i] = new TMatrixF(43,191); +// cf->WriteFullMatrix(fmatrix, i); + // cf->WriteFullNMatrix(fnmat, i); +// delete fmatrix; + // delete fnmat; //fnmat[i] = new TMatrixI(8213,8213); } else { - 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) - fnqplane[i] = new TMatrixD(43,191);// + 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; + // }; }; }; @@ -365,6 +394,8 @@ // // FULL CALORIMETER // + fqplane = cf->LoadFullAverage(rbi); + fnqplane = cf->LoadFullNAverage(rbi); CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); // Int_t nplane = 0; @@ -388,7 +419,8 @@ // 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[rbi])[nplane][mstrip] += 1.; + if ( mstrip < (191-cs) ) (*fnqplane)[nplane][mstrip] += 1.; }; }; }; @@ -407,10 +439,17 @@ // mstrip = cd + strip; // - (*fqplane[rbi])[nplane][mstrip] += mip; + // (*fqplane[rbi])[nplane][mstrip] += mip; + (*fqplane)[nplane][mstrip] += mip; // }; - + // + cf->WriteFullMean(fqplane, rbi); + cf->WriteFullNMean(fnqplane, rbi); + cf->UnLoadFullAverage(rbi); + cf->UnLoadFullNAverage(rbi); + delete fqplane; + delete fnqplane; // }; } @@ -435,26 +474,42 @@ // }; } else { - + // for (Int_t i=0; iLoadFullAverage(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]; +// 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]; } else { - (*fqplane[i])[j][k] = 0.; + (*fqplane)[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]); + 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; }; - for (Int_t i=0; iWriteFullMean(fqplane[i], i); - // - }; + // +// for (Int_t i=0; iWriteFullMean(fqplane[i], i); +// // +// }; }; // cf->WriteNumBin(nbin); @@ -485,6 +540,7 @@ }; }; // +if ( rbi != 3 ) return; if ( erig < rig[0] ) return; if ( erig >= rig[nbin-1] ) return; // @@ -527,6 +583,13 @@ // // FULL CALORIMETER // + printf(" Retrieve matrix %i IEV %i \n",rbi,iev); + // fmatrix = cf->LoadFullMatrix(rbi); +// cf->LoadFullMatrix(rbi,fmatrix); + // fnmat = cf->LoadFullNMatrix(rbi); + printf(" done \n"); + printf(" start loop \n"); + // CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); // Int_t nplane = 0; @@ -535,8 +598,14 @@ Int_t strip = 0; Float_t mip = 0.; // - Int_t mindgf = 48; - Int_t dgf = 143; + // 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 cs = 0; Int_t cd = 0; Int_t mstrip = 0; @@ -561,10 +630,10 @@ // }; // - Float_t mip1; + Float_t mip1 = 1.; Int_t cs1; Int_t cd1; - Float_t mip2; + Float_t mip2 = 1.; Int_t cs2; Int_t cd2; Int_t mi = -1; @@ -580,49 +649,82 @@ Int_t mstrip2min = 0; Int_t mstrip2max = 0; // + Int_t toto = 0; + // for (Int_t nplane1=0; nplane1<43; nplane1++){ - for (Int_t mstrip1=0; mstrip1<191; mstrip1++){ + if ( nplane1 >= 37 ) nn1 = nplane1 + 1; + vi1 = 1; + if ( nn1%2 ) vi1 = 0; + pl1 = (nn1 - 1 + vi1)/2; + // + cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1 + // + cd1 = 95 - cs1; + // + mstrip1min = TMath::Max(mindgf,(cd1+0)); + mstrip1max = TMath::Min(dgf,(cd1+95)) + 1; + // + if ( nplane1 == 0 || nplane1 == 42 ) printf(" pl %i mstrip1min %i mstrip1max %i \n",nplane1,mstrip1min,mstrip1max); + // + for (Int_t mstrip1=mstrip1min; mstrip1= 37 ) nn1 = nplane1 + 1; - vi1 = 1; - if ( nn1%2 ) vi1 = 0; - pl1 = (nn1 - 1 + vi1)/2; - // - cs1 = ct->tibar[pl1][vi1] - 1; // convertire nplane in pl1 e vi1 - // - cd1 = 95 - cs1; - // - mstrip1min = cd1 + 0; - mstrip1max = cd1 + 95; - mip1 = mipv[nplane1][mstrip1]; - // - if ( mstrip1 >= mindgf && mstrip1 <= dgf ){ - mi++; - // + mip1 = mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi); + // + mi = (nplane1 * 191) + mstrip1; + // + // if ( mstrip1 > mstrip1min ) break; + // if ( mstrip1 > dgf ) break; + // if ( mstrip1 >= mindgf && mstrip1 <= dgf && mstrip1 >= mstrip1min && mstrip1 <= mstrip1max ){ + // + // finmat[nplane1][mstrip1]++; + (*fnmat[rbi])[nplane1][mstrip1] += 1.; + // + if ( mip1 != 0. ){ + // for (Int_t nplane2=0; nplane2<43; nplane2++){ - for (Int_t mstrip2=0; mstrip2<191; mstrip2++){ - // - if ( nplane2 >= 37 ) nn2 = nplane2 + 1; - vi2 = 1; - if ( nn2%2 ) vi2 = 0; - pl1 = (nn2 - 1 + vi2)/2; - // - cs2 = ct->tibar[pl2][vi2] - 1; + // + if ( nplane2 >= 37 ) nn2 = nplane2 + 1; + vi2 = 1; + if ( nn2%2 ) vi2 = 0; + pl1 = (nn2 - 1 + vi2)/2; + // + cs2 = ct->tibar[pl2][vi2] - 1; + // + cd2 = 95 - cs2; + // + // mstrip2min = cd2 + 0; + // mstrip2max = cd2 + 95; + mstrip2min = TMath::Max(mindgf,(cd2+0)); + mstrip2max = TMath::Min(dgf,(cd2+95)) + 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); // - mstrip2min = cd2 + 0; - mstrip2max = cd2 + 95; - // - mip2 = mipv[nplane2][mstrip2]; - // - if ( mstrip2 >= mindgf && mstrip2 <= dgf ){ - mj++; - (*fmatrix[rbi])[mi][mj] += (mipv[nplane1][mstrip1] - cf->GetFullAverageAt(nplane1,mstrip1,erig)) * (mipv[nplane2][mstrip2] - cf->GetFullAverageAt(nplane2,mstrip2,erig)); - if ( mstrip1 >= mstrip1min && mstrip1 <= mstrip1max && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max){ - (*fnmat[rbi])[mi][mj] += 1.; - }; + if ( mip2 != 0. ){ + // + mj = (nplane2 * 191) + mstrip2; + // mj++; + // + // if ( mstrip2 > mstrip2min ) break; + // if ( mstrip2 > dgf ) break; + // if ( mstrip2 >= mindgf && mstrip2 <= dgf && mstrip2 >= mstrip2min && mstrip2 <= mstrip2max ){ + // 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)[mi][mj] += (mip1 * mip2) * 1000000.; + toto++; + // (*fmatrix)[mi][mj] += 1.; + // cf->GetFullAverageAt(nplane1,mstrip1,erig,rbi); + // cf->GetFullAverageAt(nplane2,mstrip2,erig,rbi); + // (*fnmat)[mi][mj] += 1.; + // }; }; }; }; @@ -630,6 +732,20 @@ }; }; // + printf(" toto = %i \n",toto); + printf("\n done \n"); + printf(" write matrix \n"); +// cf->WriteFullMatrix(fmatrix, rbi); + // cf->WriteFullNMatrix(fnmat, rbi); + printf(" done \n"); + printf(" unload matrix \n"); + // cf->UnLoadFullMatrix(rbi); + // cf->UnLoadFullNMatrix(rbi); + printf(" done \n"); + printf(" delete matrix \n"); +// delete fmatrix; + // delete fnmat; + printf(" done \n"); }; } @@ -685,30 +801,142 @@ // // FULL // - for (Int_t i=0; i 0. ){ - (*fmatrix[i])[ii][j] /= (*fnmat[i])[ii][j]; - } else { - (*fmatrix[i])[ii][j] = 0.; +// 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.; +// }; +// }; +// }; +// +// 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; + // j1 = -1; + for (Int_t iij=0; iij<43; iij++){ + for (Int_t jj=0; jj<191; jj++){ + // + j1 = (iij * 191) + 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] = 1.; + } else { + (*fmatrix)[i1][j1] /= (*fnmat[i])[ii][j]; + nonzero++; + if ( i1 == 0 ) nonzero1++; + }; + // 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]; + }; + + }; }; }; }; // - cf->WriteFullMatrix(fmatrix[i],i); + printf(" Matrix has %i non-zero elements \n",nonzero); + 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"); + // + // + // 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[i]->Determinant() == 0. ){ + 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)); + // 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 ( mymat->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)(mymat->Invert(&det)); + printf(" Bin %i determinant is %f \n",i,det); + cf->WriteInvertedFullMatrix((TMatrixF)invmatrix,i); + }; + cf->WriteFullMatrix(mymat, 99); + + + // + cf->UnLoadFullMatrix(i); + // cf->UnLoadFullNMatrix(i); + delete fmatrix; + // delete fnmat; + // }; }; //