--- calo/flight/CaloFranzini/src/Calib.cpp 2007/12/13 17:08:09 1.1 +++ calo/flight/CaloFranzini/src/Calib.cpp 2007/12/18 09:55:05 1.2 @@ -14,17 +14,18 @@ // using namespace std; // +extern Bool_t MATRIX; CaloFranzini *cf; Int_t nbin; Float_t rig[18]; -Float_t rmean[18]; -Float_t qqplane[18][43]; -Int_t nnqplane[18][43]; - -TArrayF *qplane[18]; -TArrayI *nqplane[18]; -TMatrixD *matrix[18]; -TMatrixD *nmat[18]; +Float_t rmean[17]; +//Float_t qqplane[17][43]; +//Int_t nnqplane[17][43]; + +TArrayF *qplane[17]; +TArrayI *nqplane[17]; +TMatrixD *matrix[17]; +TMatrixD *nmat[17]; //=============================================================================== bool Select( PamLevel2* event ){ @@ -134,7 +135,14 @@ void CreateHistos( PamLevel2* event , TString file){ cf = new CaloFranzini(event); - cf->CreateMatrixFile(file.Data()); + // + if ( MATRIX ){ + cf->UpdateMatrixFile(file.Data()); + cf->LoadBin(); + } else { + cf->CreateMatrixFile(file.Data()); + }; + // // nbin = 18; rig[0] = 0.1; @@ -156,47 +164,30 @@ rig[16] = 200.; rig[17] = 400.; // -// memset(qplane, 0, 44*18*sizeof(Float_t)); - memset(rmean, 0, 18*sizeof(Float_t)); - memset(qqplane, 0, 43*18*sizeof(Float_t)); - memset(nnqplane, 0, 43*18*sizeof(Float_t)); + memset(rmean, 0, 17*sizeof(Float_t)); // - for (Int_t i=0; i < 18 ; i++){ + for (Int_t i=0; i < 17 ; i++){ matrix[i] = new TMatrixD(43,43); qplane[i] = new TArrayF(43); nqplane[i] = new TArrayI(43); nmat[i] = new TMatrixD(43,43); }; -// for (Int_t i=0; i < 18 ; i++){ -// for (Int_t j = 0; j < 43;i++){ -// (*qplane[i])[j] = 0.; -// qplane[i]->Reset(); -// nqplane[i]->Reset(); -// for (Int_t k = 0; k < 43;k++){ -// (*matrix[i])[j][k] = 0.; -// (*nmat[i])[j][k] = 0.; -// }; -// }; -// }; // } //=============================================================== void FindAverage( PamLevel2* L2, int iev ){ // - // printf("SELEZIONATO \n"); Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity(); // Int_t rbi = 0; - for (Int_t i = 0; i=rig[i] && erig < rig[i+1] ){ rbi = i; break; }; }; // - // printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]); - // if ( erig < rig[0] ) return; if ( erig >= rig[nbin-1] ) return; // @@ -204,27 +195,10 @@ // Int_t dgf = 43; // -// for (Int_t i=0; i < 22; i++){ -// if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){ -// dgf = 2 * i; - // break; - // }; - // if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){ - // dgf = 1 + 2 * i; - // break; - // }; - // }; - // -// if ( dgf < 43 && dgf > 37 ) dgf--; - // for (Int_t i=0; i 37 ) nplane--; - // printf(" mip %f view %i plane %i => nplane %i \n",mip,view,plane,nplane); if ( nplane < dgf ){ - // printf(" PRIMA: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]); (*qplane[rbi])[nplane] += mip; - qqplane[rbi][nplane] += mip; - // printf(" DOPO: qpl %f nqpl %i \n",(*qplane[rbi])[nplane],(*nqplane[rbi])[nplane]); }; // }; } +void CalculateAverage(){ + for (Int_t i=0; i 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0]; + for (Int_t j=0; j<43 ; j++){ + if ( (*nqplane[i])[j] > 0 ){ + (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j]; + } else { + (*qplane[i])[j] = 0.; + }; + printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]); + }; + }; + + cf->WriteNumBin(nbin); + + TArrayF *rigbin = new TArrayF(18, rig); + cf->WriteRigBin(rigbin); + + TArrayF *rmeanbin = new TArrayF(17, rmean); + TFile *file = cf->GetFile(); + file->cd(); + file->WriteObject(&(*rmeanbin), "binrigmean"); + // + for (Int_t i=0; iWriteLongMean(qplane[i], i); + // + }; + // +} + //=============================================================== void FindMatrix( PamLevel2* L2, int iev ){ // - // printf("SELEZIONATO \n"); Float_t erig = L2->GetTrack(0)->GetTrkTrack()->GetRigidity(); // Int_t rbi = 0; @@ -264,34 +264,18 @@ }; }; // - // printf("SELEZIONATO rbi %i erig %f rig[0] %f rig[nbin] %f \n",rbi, erig, rig[0],rig[nbin-1]); - // if ( erig < rig[0] ) return; if ( erig >= rig[nbin-1] ) return; // Int_t dgf = 43; - // -// for (Int_t i=0; i < 22; i++){ -// if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][1] < 0 ){ -// dgf = 2 * i; - // break; - // }; - // if ( L2->GetTrack(0)->GetCaloTrack()->tibar[i][0] < 0 ){ - // dgf = 1 + 2 * i; - // break; - // }; -// }; - // -// if ( dgf < 43 && dgf > 37 ) dgf--; +// printf("1qui \n"); // for (Int_t i=0; iGetAverageAt(i,erig)) * (hpl[j] - cf->GetAverageAt(j,erig)); }; }; +// printf("4qui \n"); + // gObjectTable->Print(); } -void CalculateAverage(){ - for (Int_t i=0; i 0 ) rmean[i] /= (Float_t)(*nqplane[i])[0]; - for (Int_t j=0; j<43 ; j++){ - if ( (*nqplane[i])[j] > 0 ){ - (*qplane[i])[j] /= (Float_t)(*nqplane[i])[j]; - } else { - (*qplane[i])[j] = 0.; - }; - printf(" BIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],(*qplane[i])[j],(*nqplane[i])[j]); -// if ( nnqplane[i][j] > 0 ){ -// qqplane[i][j] /= (Float_t)nnqplane[i][j]; -// } else { -// qqplane[i][j] = 0.; -// }; -// printf(" BBIN %i plane %i average energy %f qplane %f nqplane %i \n",i,j,rmean[i],qqplane[i][j],nnqplane[i][j]); - }; - }; - - cf->WriteNumBin(nbin); - - TArrayF *rigbin = new TArrayF(18, rig); - cf->WriteRigBin(rigbin); - - TArrayF *rmeanbin = new TArrayF(18, rmean); - TFile *file = cf->GetFile(); - file->cd(); - // rigbin->Write("binrig"); - file->WriteObject(&(*rmeanbin), "binrigmean"); - -// cf->WriteRigBin(rigbin); - - for (Int_t i=0; iWriteLongMean(qplane[i], i); - - }; -} //=============================================================== // Save histograms // @@ -370,70 +318,51 @@ // //=============================================================== void SaveHistos(){ - - printf("Finished, calculating average and inverting matrices\n"); - - - -// cf->WriteNumBin(nbin); - -// TArrayF *rigbin = new TArrayF(18, rig); -// cf->WriteRigBin(rigbin); - - for (Int_t i=0; iWriteLongMean(qplane[i], i); - + // + if ( MATRIX ){ + // + printf("Finished, calculating average and inverting matrices\n"); // - // determine the average matrix - // - for (Int_t ii=0; ii<43; ii++){ - for (Int_t j=0; j<43; j++){ - // if ( i == 1 ) printf(" num %f \n",(*nmat[i])[ii][j]); - if ( (*nmat[i])[ii][j] > 0. ){ - // if ( i == 1 ) printf(" value %f \n",(*matrix[i])[ii][j]); - (*matrix[i])[ii][j] /= (*nmat[i])[ii][j]; - } else { - (*matrix[i])[ii][j] = 0.; + for (Int_t i=0; i 0. ){ + (*matrix[i])[ii][j] /= (*nmat[i])[ii][j]; + } else { + (*matrix[i])[ii][j] = 0.; + }; }; }; - }; - -// if ( i == 2 ){ -// printf("\n"); -// for (Int_t ii=0; ii<43; ii++){ -// for (Int_t j=0; j<43; j++){ -// printf(" %.f",(*matrix[i])[ii][j]); -// }; -// printf("\n"); -// }; -// printf("\n"); -// }; - - - cf->WriteLongMatrix(matrix[i],i); - - if ( matrix[i]->Determinant() == 0. ){ - printf("\n"); - for (Int_t ii=0; ii<43; ii++){ - for (Int_t j=0; j<43; j++){ - printf(" %.f",(*matrix[i])[ii][j]); - }; - printf("\n"); + // + cf->WriteLongMatrix(matrix[i],i); + // + if ( matrix[i]->Determinant() == 0. ){ + printf("\n"); + for (Int_t ii=0; ii<43; ii++){ + for (Int_t j=0; j<43; j++){ + printf(" %.f",(*matrix[i])[ii][j]); + }; + printf("\n"); + }; + printf("\n"); + 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)(matrix[i]->Invert(&det)); + printf(" Bin %i determinant is %f \n",i,det); + cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i); }; - printf("\n"); - 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)(matrix[i]->Invert(&det)); - printf(" Bin %i determinant is %f \n",i,det); - cf->WriteInvertedLongMatrix((TMatrixD)invmatrix,i); }; + // + printf(" done, closing file and exiting\n"); + // }; - - printf(" done, closing file and exiting\n"); - + // cf->CloseMatrixFile(); - + // cf->Delete(); + // }