--- calo/flight/CaloFranzini/src/CaloFranzini.cpp 2007/12/18 09:55:07 1.4 +++ calo/flight/CaloFranzini/src/CaloFranzini.cpp 2009/12/14 14:35:53 1.11 @@ -16,7 +16,8 @@ CaloFranzini::CaloFranzini(PamLevel2 *l2p){ // - file = NULL; + lfile = NULL; + ffile = NULL; brig = NULL; brigm = NULL; nbin = 0; @@ -31,7 +32,13 @@ dolong = true; dofull = false; sntr = 0; + OBT = 0; + PKT = 0; + atime = 0; // + usepl18x = false; + // + crig = false; sel = true; cont = false; N = 0; @@ -45,15 +52,22 @@ void CaloFranzini::Clear(){ // - OBT = 0; - PKT = 0; - atime = 0; - longtzeta = -100.; - fulltzeta = -100.; + longtzeta = 0.; + fulltzeta = 0.; degfre = 0; + fdegfre = 0; memset(estrip, 0, 2*22*96*sizeof(Float_t)); memset(qplane, 0, 43*sizeof(Float_t)); // + numneg = 0; + numpos = 0; + negfulltzeta = 0.; + posfulltzeta = 0.; + aveposvar = 0.; + avenegvar = 0.; + minsvalue = numeric_limits::max(); + maxsvalue = numeric_limits::min(); + // } void CaloFranzini::Print(){ @@ -63,25 +77,37 @@ printf("========================================================================\n"); printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); printf(" debug [debug flag]:.. %i\n",debug); - printf(" degree of freedom :.. %i\n",this->GetDegreeOfFreedom()); + printf(" long degree of freedom :.. %i\n",this->GetLongDegreeOfFreedom()); printf(" longtzeta :.. %f\n",longtzeta); printf(" longtzeta normalized :.. %f\n",this->GetNormLongTZeta()); - // printf(" fulltzeta :.. %f\n",fulltzeta); - // printf(" longtzeta normalized :.. %f\n",this->GetNormFullTZeta()); + printf(" full degree of freedom :.. %i\n",this->GetFullDegreeOfFreedom()); + printf(" fulltzeta :.. %f\n",fulltzeta); + printf(" fulltzeta normalized :.. %f\n",this->GetNormFullTZeta()); + printf(" fulltzeta negative contribute :.. %f\n",negfulltzeta); + printf(" fulltzeta positive contribute :.. %f\n",posfulltzeta); + printf(" fulltzeta number of negatives :.. %i\n",numneg); + printf(" fulltzeta number of positives :.. %i\n",numpos); + printf(" fulltzeta minimum variance :.. %f\n",minsvalue); + printf(" fulltzeta maximum variance :.. %f\n",maxsvalue); + printf(" fulltzeta average positive var :.. %f\n",aveposvar); + printf(" fulltzeta average negative var :.. %f\n",avenegvar); printf("========================================================================\n"); // } void CaloFranzini::Delete(){ // - if ( file ) file->Close(); + if ( ffile ) ffile->Close(); + if ( lfile ) lfile->Close(); // Clear(); // } void CaloFranzini::SetNoWpreSampler(Int_t n){ - if ( NC+n < 23 ){ + Int_t nc2 = NC/2; + if ( NC >= 37 ) nc2 = (NC+1)/2; + if ( nc2+n < 23 ){ N = n; } else { printf(" ERROR! Calorimeter is made of 22 W planes\n"); @@ -105,6 +131,19 @@ }; } +void CaloFranzini::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){ + this->SetNoWpreSampler(0); + this->SetNoWcalo(0); + if ( NoWpreSampler < NoWcalo ){ + this->SetNoWpreSampler(NoWpreSampler); + this->SetNoWcalo(NoWcalo); + } else { + this->SetNoWcalo(NoWcalo); + this->SetNoWpreSampler(NoWpreSampler); + }; +} + + void CaloFranzini::Process(){ this->Process(0); } @@ -117,7 +156,7 @@ return; }; // - if ( !nbin || !file || !brig ){ + if ( !nbin || !brig || (!lfile && !ffile) ){ printf(" ERROR: it seems covariance matrix file has not been opened (CaloFranzini::Open()) \n"); printf(" ERROR: CaloFranzini variables _NOT_ filled \n"); return; @@ -146,6 +185,7 @@ this->Clear(); // Float_t rig = L2->GetTrack(itr)->GetTrkTrack()->GetRigidity(); + if ( crig ) rig = L2->GetCaloLevel2()->qtot/260.; // // Fill the estrip matrix // @@ -154,24 +194,27 @@ Int_t plane = 0; Int_t strip = 0; Float_t mip = 0.; - for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ - // - mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); - // - estrip[view][plane][strip] = mip; - // - nplane = 1 - view + 2 * (plane - N); - if ( nplane > (37-(2*N)) ) nplane--; - // - if ( plane == (18+N) ) mip = 0.; - if ( nplane > -1 ) qplane[nplane] += mip; - // - }; // // // if ( dolong ){ // + for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + if ( !usepl18x && view==0 && plane==18 ) mip = 0.; + // + estrip[view][plane][strip] = mip; + // + nplane = 1 - view + 2 * (plane - N); + if ( nplane > (37-(2*N)) ) nplane--; + // + // if ( plane == (18+N) ) mip = 0.; + if ( nplane > -1 ) qplane[nplane] += mip; + // + }; + // if ( cont ){ for (Int_t i=0; i<22; i++){ if ( i == (18+N) ){ @@ -209,24 +252,250 @@ // degfre = TMath::Min(dgf,NC); // + // Float_t longzdiag = 0.; + // Float_t longzout = 0.; + // if ( degfre > 0 ){ for (Int_t i = 0; i < degfre; i++){ if ( i != mask18b ){ for (Int_t j = 0; j < degfre; j++){ if ( j != mask18b ){ +// if ( i == j ){ +// longzdiag += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig)); +// if ( debug ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig)); +// } else { +// longzout += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig)); +// if ( debug && i == (j+1) ) printf(" %i %i %f %f %f %f %f\n",i,j,qplane[i],this->GetAverageAt(i,rig),this->GetHmatrixAt(i,j,rig),qplane[j],this->GetAverageAt(j,rig)); +// }; + // longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig)); + // }; }; }; }; + // if ( debug ) printf(" diagonal total %f out of diagonal total %f \n",longzdiag,longzout); }; }; if ( dofull ){ - printf(" ERROR: NOT IMPLEMENTED YET\n"); - printf(" ERROR: CaloFranzini variables _NOT_ filled \n"); - }; + // printf(" ERROR: NOT IMPLEMENTED YET\n"); + // printf(" ERROR: CaloFranzini variables _NOT_ filled \n"); + // + // FULL CALORIMETER + // + CaloTrkVar *ct = L2->GetTrack(0)->GetCaloTrack(); + // + Int_t dgf = 0; + Int_t cs = 0; + Int_t cd = 0; + Int_t mstrip = 0; + // + Int_t maxnpl = 0; + Float_t mipv[43][11]; + memset(mipv,0,43*11*sizeof(Float_t)); + // + for ( Int_t i=0; iGetCaloLevel1()->istrip; i++ ){ + // + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + // + if ( !usepl18x && view==0 && plane==18 ) mip = 0.; + // + // nplane = 1 - view + 2 * plane; + // if ( nplane > 37 ) nplane--; + nplane = 1 - view + 2 * (plane - N); + if ( nplane > (37-(2*N)) ) nplane--; + // + cs = ct->tibar[plane][view] - 1; + // + if ( ct->tibar[plane][view] > -1 ){ + // + cd = 95 - cs; + // + mstrip = cd + strip; + // + Int_t lstr = this->ConvertStrip(mstrip); + // + if ( nplane > maxnpl ) maxnpl = nplane; + if ( nplane > -1 ) mipv[nplane][lstr] += mip; + // + }; + // mipv[nplane][lstr] += mip; + // + }; + // + Float_t mip1 = 1.; + Int_t cs1; + Int_t cd1; + Float_t mip2 = 1.; + Int_t cs2; + Int_t cd2; + Int_t mi = -1; + Int_t mj = -1; + Int_t nn1 = 0; + Int_t pl1 = 0; + Int_t vi1 = 0; + Int_t nn2 = 0; + Int_t pl2 = 0; + Int_t vi2 = 0; + Int_t mstrip1min = 0; + Int_t mstrip1max = 0; + Int_t mstrip2min = 0; + Int_t mstrip2max = 0; + // + Int_t tonpl = TMath::Min(maxnpl,NC); + // + Int_t rbi = 0; + for (Int_t i = 0; i=brigm->At(i) && rig < brigm->At(i+1) ){ + rbi = i; + break; + }; + }; + // + // + Int_t therigb = rbi; + // + if ( rig < brigm->At(2) ){ +// therigb = 0; + therigb = 2; + }; + // if ( rig >= brigm->At(nbin-5) ){ + if ( rig >= brigm->At(nbin-2) ){ + therigb = nbin - 3; + // therigb = nbin - 5; + }; + Int_t mtherig = therigb; + if ( mtherig >= 13 ) mtherig = 12; + // + if ( debug ) printf(" rig %f brigm0 %f brigm n2 %f nbin %i \n",rig,brigm->At(0),brigm->At(nbin-2),nbin); + // + if ( cont ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18+N) ){ + mask18b = 1 + 2 * (i - N); + break; + }; + }; + }; + // + if ( sel ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18-N) ){ + mask18b = 1 + 2 * (i - N); + break; + }; + }; + }; + // + if ( mask18b == 37 ) mask18b = -1; + // + if ( debug ) printf("Start main loop \n"); + // + for (Int_t nplane1=0; nplane1= 37 ) nn1 = nplane1 + 1; + vi1 = 1; + if ( nn1%2 ) vi1 = 0; + pl1 = (nn1 - 1 + vi1)/2; + // + cs1 = ct->tibar[pl1][vi1] - 1; + // + if ( ct->tibar[pl1][vi1] > -1 ){ + // + dgf++; + // + cd1 = 95 - cs1; + // + Int_t at1 = TMath::Max(0,(cd1+0)); + Int_t at2 = TMath::Min(190,(cd1+95)); + mstrip1min = this->ConvertStrip(at1); + mstrip1max = this->ConvertStrip(at2) + 1; + // + for (Int_t mstrip1=mstrip1min; mstrip1GetFullAverageAt(nplane1,mstrip1,rig,therigb); + if ( mip1 < 0 ){ + numneg++; + avenegvar += mip1; + }; + if ( mip1 >= 0 ){ + numpos++; + aveposvar += mip1; + }; + if ( mip1 < minsvalue ) minsvalue = mip1; + if ( mip1 >= maxsvalue ) maxsvalue = mip1; + // + mi = (nplane1 * 11) + mstrip1; + // + if ( mip1 != 0. ){ + // + for (Int_t nplane2=0; nplane2= 37 ) nn2 = nplane2 + 1; + vi2 = 1; + if ( nn2%2 ) vi2 = 0; + pl1 = (nn2 - 1 + vi2)/2; + // + cs2 = ct->tibar[pl2][vi2] - 1; + // + if ( ct->tibar[pl2][vi2] > -1 ){ + // + cd2 = 95 - cs2; + // + Int_t t1 = TMath::Max(0,(cd2+0)); + Int_t t2 = TMath::Min(190,(cd2+95)); + mstrip2min = this->ConvertStrip(t1); + mstrip2max = this->ConvertStrip(t2) + 1; + // + for (Int_t mstrip2=mstrip2min; mstrip2GetFullAverageAt(nplane2,mstrip2,rig,therigb); + // + if ( mip2 != 0. ){ + // +// 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; + // + Float_t svalue = mip1 * this->GetFullHmatrixAt(mi,mj,rig,mtherig,therigb) * mip2; + // fulltzeta += mip1 * this->GetFullHmatrixAt(mi,mj,rig,therigb) * mip2; + fulltzeta += svalue; + if ( svalue < 0. ) negfulltzeta += svalue; + if ( svalue > 0. ) posfulltzeta += svalue; + // (*fmatrix[rbi])[mi][mj] += (mip1 * mip2); + // if ( mstrip1 == mstrip2 ) printf("\n\n=> nplane1 %i nplane2 %i mstrip1 %i mstrip2 %i \n => mip1 %f H %f mip2 %f \n => mipv1 %f ave1 %f mipv2 %f ave2 %f \n => rig %f rigbin %i mtherigb %i\n",nplane1,nplane2,mstrip1,mstrip2,mip1,this->GetFullHmatrixAt(mi,mj,rig,mtherig),mip2,mipv[nplane1][mstrip1],this->GetFullAverageAt(nplane1,mstrip1,rig,therigb),mipv[nplane2][mstrip2],this->GetFullAverageAt(nplane2,mstrip2,rig,therigb),rig,therigb,mtherig); + // + }; + }; + }; + }; + }; + }; + }; + }; + }; + }; + // + fdegfre = dgf*11; + if ( numpos ) aveposvar /= numpos; + if ( numneg ) avenegvar /= numneg; + // + }; // - if ( debug ) this->Print(); + // if ( debug ) this->Print(); if ( debug ) printf(" exit \n"); } @@ -234,25 +503,32 @@ Float_t CaloFranzini::GetNormLongTZeta(){ Process(); Float_t normz = 0.; - if ( degfre > 0 ) normz = longtzeta/(Float_t)degfre; + if ( degfre != 0 ) normz = longtzeta/(Float_t)degfre; return normz; } Float_t CaloFranzini::GetNormFullTZeta(){ Process(); Float_t normz = 0.; - if ( degfre > 0 ) normz = fulltzeta/(Float_t)degfre; + if ( fdegfre != 0 ) normz = fulltzeta/(Float_t)fdegfre; return normz; } Bool_t CaloFranzini::CreateMatrixFile(TString matrixfile){ // - file = new TFile(matrixfile.Data(),"READ"); + lfile = new TFile(matrixfile.Data(),"READ"); // - if ( !file || file->IsZombie() ){ - file = new TFile(matrixfile.Data(),"RECREATE"); - printf(" Create file %s \n",file->GetName()); + if ( !lfile || lfile->IsZombie() ){ + if ( dolong ){ + lfile = new TFile(matrixfile.Data(),"RECREATE"); + printf(" Create file %s \n",lfile->GetName()); + }; + if ( dofull ){ + ffile = new TFile(matrixfile.Data(),"RECREATE"); + printf(" Create file %s \n",ffile->GetName()); + }; } else { + lfile->Close(); printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data()); return(false); }; @@ -263,11 +539,22 @@ Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){ // - file = new TFile(matrixfile.Data(),"UPDATE"); + if ( dolong ){ + lfile = new TFile(matrixfile.Data(),"UPDATE"); + // + if ( !lfile || lfile->IsZombie() ){ + printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data()); + return(false); + }; + }; // - if ( !file || file->IsZombie() ){ - printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data()); - return(false); + if ( dofull ){ + ffile = new TFile(matrixfile.Data(),"UPDATE"); + // + if ( !ffile || ffile->IsZombie() ){ + printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data()); + return(false); + }; }; // return(true); @@ -275,66 +562,165 @@ } void CaloFranzini::WriteNumBin(Int_t numbin){ - file->cd(); - TArrayI nbi(1, &numbin); - file->WriteObject(&nbi, "nbinenergy"); + if ( dolong ){ + lfile->cd(); + TArrayI nbi(1, &numbin); + lfile->WriteObject(&nbi, "nbinenergy"); + }; + if ( dofull ){ + ffile->cd(); + TArrayI nbi(1, &numbin); + ffile->WriteObject(&nbi, "nbinenergy"); + }; } void CaloFranzini::WriteRigBin(TArrayF *rigbin){ - file->cd(); - // rigbin->Write("binrig"); - file->WriteObject(&(*rigbin), "binrig"); + if ( dolong ){ + lfile->cd(); + // rigbin->Write("binrig"); + lfile->WriteObject(&(*rigbin), "binrig"); + }; + if ( dofull ){ + ffile->cd(); + // rigbin->Write("binrig"); + ffile->WriteObject(&(*rigbin), "binrig"); + }; } void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){ - file->cd(); + lfile->cd(); TString name = Form("qplmeann%i",bin); - file->WriteObject(&(*qpl),name.Data()); + lfile->WriteObject(&(*qpl),name.Data()); +} + +void CaloFranzini::WriteFullMean(TMatrixD *qpl, Int_t bin){ + ffile->cd(); + TString name = Form("fqplmeann%i",bin); + ffile->WriteObject(&(*qpl),name.Data()); + ffile->Purge(); +} + +void CaloFranzini::WriteFullNMean(TMatrixD *qpl, Int_t bin){ + ffile->cd(); + TString name = Form("fnqplmeann%i",bin); + ffile->WriteObject(&(*qpl),name.Data()); + ffile->Purge(); } void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){ - file->cd(); + lfile->cd(); TString name = Form("matrixn%i",bin); // mat.Write(name.Data()); - file->WriteObject(&mat,name.Data()); + lfile->WriteObject(&mat,name.Data()); +} + +void CaloFranzini::WriteInvertedFullMatrix(TMatrixD mat, Int_t bin){ + ffile->cd(); + TString name = Form("fmatrixn%i",bin); + // mat.Write(name.Data()); + ffile->WriteObject(&mat,name.Data()); } void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){ - file->cd(); + lfile->cd(); TString name = Form("origmatrixn%i",bin); // mat.Write(name.Data()); - file->WriteObject(&(*mat),name.Data()); + lfile->WriteObject(&(*mat),name.Data()); } -void CaloFranzini::CloseMatrixFile(){ - file->cd(); - file->Close(); +void CaloFranzini::WriteFullMatrix(TMatrixD *mat, Int_t bin){ + ffile->cd(); + TString name = Form("origfmatrixn%i",bin); + // mat.Write(name.Data()); + ffile->WriteObject(&(*mat),name.Data()); + ffile->Purge(); } +void CaloFranzini::WriteFullNMatrix(TMatrixF *mat, Int_t bin){ + ffile->cd(); + TString name = Form("origfnmatrixn%i",bin); + // mat.Write(name.Data()); + ffile->WriteObject(&(*mat),name.Data()); + ffile->Purge(); +} + +void CaloFranzini::CloseMatrixFile(){ + if ( dolong ){ + lfile->cd(); + lfile->Close(); + }; + if ( dofull ){ + ffile->cd(); + ffile->Close(); + }; +} Bool_t CaloFranzini::Open(TString matrixfile){ + return(this->Open(matrixfile,"")); +} + +Bool_t CaloFranzini::Open(TString longmatrixfile, TString fullmatrixfile){ // // find matrix file // - if ( !strcmp(matrixfile.Data(),"") ){ - if (dolong) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root"; - if (dofull) matrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root"; + if ( !strcmp(longmatrixfile.Data(),"") ){ + if (dolong) longmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_longel.root"; }; - // - file = new TFile(matrixfile.Data(),"READ"); - // - if ( !file || file->IsZombie() ){ - printf(" ERROR: cannot open file %s \n",matrixfile.Data()); - return(false); + if ( !strcmp(fullmatrixfile.Data(),"") ){ + if (dofull) fullmatrixfile = (TString)gSystem->ExpandPathName("$PAM_CALIB")+"/cal-param/covmatrix_fullel.root"; }; // - if ( !this->LoadBin() ){ - printf(" %s \n",matrixfile.Data()); - return(false); + if ( dolong ){ + // + lfile = new TFile(longmatrixfile.Data(),"READ"); + // + if ( !lfile || lfile->IsZombie() ){ + printf(" ERROR: cannot open file %s \n",longmatrixfile.Data()); + return(false); + }; + // + if ( !this->LoadBin() ){ + printf(" %s \n",longmatrixfile.Data()); + return(false); + }; + // + if ( !this->LoadLong() ){ + printf(" %s \n",longmatrixfile.Data()); + return(false); + }; + // + if ( !this->LoadMatrices() ){ + printf(" %s \n",longmatrixfile.Data()); + return(false); + }; }; - if ( !this->LoadMatrices() ){ - printf(" %s \n",matrixfile.Data()); - return(false); + // + if ( dofull ){ + // + ffile = new TFile(fullmatrixfile.Data(),"READ"); + // + if ( !ffile || ffile->IsZombie() ){ + printf(" ERROR: cannot open file %s \n",fullmatrixfile.Data()); + return(false); + }; + // + if ( !dolong ){ + // + if ( !this->LoadBin(true) ){ + printf(" %s \n",fullmatrixfile.Data()); + return(false); + }; + // + }; + if ( !this->LoadFull() ){ + printf(" %s \n",fullmatrixfile.Data()); + return(false); + }; + // + if ( !this->LoadFullMatrices() ){ + printf(" %s \n",fullmatrixfile.Data()); + return(false); + }; }; // // @@ -342,38 +728,90 @@ // } + Bool_t CaloFranzini::LoadBin(){ + return(this->LoadBin(false)); +} + +Bool_t CaloFranzini::LoadBin(Bool_t full){ // - TArrayI *numbin = (TArrayI*)file->Get("nbinenergy"); - if ( !numbin ){ - printf(" ERROR: cannot read number of bins from file "); - return(false); - }; - nbin = (Int_t)numbin->At(0); - if ( nbin <= 0 ){ - printf(" ERROR: cannot work with 0 energy bins from file "); - return(false); + if ( !full ) { + // + TArrayI *numbin = (TArrayI*)lfile->Get("nbinenergy"); + if ( !numbin ){ + printf(" ERROR: cannot read number of bins from file "); + return(false); + }; + nbin = (Int_t)numbin->At(0); + if ( nbin <= 0 ){ + printf(" ERROR: cannot work with 0 energy bins from file "); + return(false); + }; + // + brig = (TArrayF*)lfile->Get("binrig"); + if ( !brig ){ + printf(" ERROR: cannot read rigidity binning from file "); + return(false); + }; + // + brigm=(TArrayF*)lfile->Get("binrigmean"); + if ( !brigm ){ + printf(" ERROR: cannot read mean rigidity binning from file "); + return(false); + }; + // + } else { + // + TArrayI *numbin = (TArrayI*)ffile->Get("nbinenergy"); + if ( !numbin ){ + printf(" ERROR: cannot read number of bins from file "); + return(false); + }; + nbin = (Int_t)numbin->At(0); + if ( nbin <= 0 ){ + printf(" ERROR: cannot work with 0 energy bins from file "); + return(false); + }; + // + brig = (TArrayF*)ffile->Get("binrig"); + if ( !brig ){ + printf(" ERROR: cannot read rigidity binning from file "); + return(false); + }; + // + brigm=(TArrayF*)ffile->Get("binrigmean"); + if ( !brigm ){ + printf(" ERROR: cannot read mean rigidity binning from file "); + return(false); + }; }; // - brig = (TArrayF*)file->Get("binrig"); - if ( !brig ){ - printf(" ERROR: cannot read rigidity binning from file "); - return(false); - }; + return(true); +} + +Bool_t CaloFranzini::LoadLong(){ // - brigm=(TArrayF*)file->Get("binrigmean"); - if ( !brigm ){ - printf(" ERROR: cannot read mean rigidity binning from file "); - return(false); + for (Int_t i=0;i<17;i++){ + TString name = Form("qplmeann%i",i); + qplmean[i] = (TArrayF*)lfile->Get(name.Data()); + if ( !qplmean[i] ){ + printf(" ERROR: cannot read average from file "); + return(false); + }; }; // + return(true); +} + +Bool_t CaloFranzini::LoadFull(){ + // for (Int_t i=0;i<17;i++){ - TString name = Form("qplmeann%i",i); - qplmean[i] = (TArrayF*)file->Get(name.Data()); - if ( !qplmean[i] ){ - printf(" ERROR: cannot read average from file "); - return(false); - }; + TString name = Form("fqplmeann%i",i); + fqplmean[i] = (TMatrixD*)ffile->Get(name.Data()); + if ( !fqplmean[i] ){ + printf(" ERROR: cannot read average from file "); + return(false); + }; }; // return(true); @@ -383,7 +821,17 @@ // for (Int_t i=0;i<17;i++){ TString name1 = Form("matrixn%i",i); - hmat[i] = (TMatrixD*)file->Get(name1.Data()); + hmat[i] = (TMatrixD*)lfile->Get(name1.Data()); + }; + // + return(true); +} + +Bool_t CaloFranzini::LoadFullMatrices(){ + // + for (Int_t i=0;i<17;i++){ + TString name1 = Form("fmatrixn%i",i); + hfmat[i] = (TMatrixD*)ffile->Get(name1.Data()); }; // return(true); @@ -411,6 +859,77 @@ // } +TMatrixD *CaloFranzini::LoadFullAverage(Int_t rigbin){ + // + TString name = Form("fqplmeann%i",rigbin); + TMatrixD *fmean=(TMatrixD*)ffile->Get(name.Data()); + // + return(fmean); + // +} + +TMatrixF *CaloFranzini::LoadFullMatrix(Int_t rigbin){ + // + TString name = Form("origfmatrixn%i",rigbin); + TMatrixF *fmatri=(TMatrixF*)ffile->Get(name.Data()); + // + return(fmatri); + // +} + +void CaloFranzini::LoadFullMatrix(Int_t rigbin, TMatrixF *fmatri){ + // + TString name = Form("origfmatrixn%i",rigbin); + fmatri=(TMatrixF*)ffile->Get(name.Data()); + // +} + +void CaloFranzini::UnLoadFullAverage(Int_t rigbin){ + // + TString name = Form("fqplmeann%i",rigbin); + ffile->Delete(name.Data()); + // +} + +void CaloFranzini::UnLoadFullMatrix(Int_t rigbin){ + // + TString name = Form("origfmatrixn%i",rigbin); + ffile->Delete(name.Data()); + // +} + +TMatrixF *CaloFranzini::LoadFullNMatrix(Int_t rigbin){ + // + TString name = Form("origfnmatrixn%i",rigbin); + TMatrixF *fnmatri=(TMatrixF*)ffile->Get(name.Data()); + // + return(fnmatri); + // +} + +void CaloFranzini::UnLoadFullNMatrix(Int_t rigbin){ + // + TString name = Form("origfnmatrixn%i",rigbin); + ffile->Delete(name.Data()); + // +} + +TMatrixD *CaloFranzini::LoadFullNAverage(Int_t rigbin){ + // + TString name = Form("fnqplmeann%i",rigbin); + TMatrixD *fnmean=(TMatrixD*)ffile->Get(name.Data()); + // + return(fnmean); + // +} + +void CaloFranzini::UnLoadFullNAverage(Int_t rigbin){ + // + TString name = Form("fnqplmeann%i",rigbin); + ffile->Delete(name.Data()); + // +} + TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){ // @@ -450,13 +969,13 @@ // if ( rig < brigm->At(0) ){ if ( rig < brig->At(0) ){ - printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0)); +// printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0)); }; therigb = 0; }; if ( rig >= brigm->At(nbin-2) ){ if ( rig >= brig->At(nbin-2) ) { - printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2)); +// printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2)); }; therigb = nbin - 3; }; @@ -474,6 +993,72 @@ Float_t b = log(maxene/minene)/(maxrig-minrig); Float_t a = minene/exp(minrig*b); Float_t ave = a*exp(b*rig); + if ( a == 0. || minene == 0. || ave != ave ){ + // if ( a == 0. || minene == 0. ){ + Float_t m = (maxene-minene)/(maxrig-minrig); + Float_t q = minene - m * minrig; + ave = rig * m + q; + }; + // + return(ave); + // +} + +Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig){ + // + Int_t therigb = 0; + for (Int_t i = 0; i=brigm->At(i) && rig < brigm->At(i+1) ){ + therigb = i; + break; + }; + }; + // + // + if ( rig < brigm->At(0) ){ + if ( rig < brig->At(0) ){ +// printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0)); + }; + therigb = 0; + }; + if ( rig >= brigm->At(nbin-2) ){ + if ( rig >= brig->At(nbin-2) ) { +// printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2)); + }; + therigb = nbin - 3; + }; + // + return(this->GetFullAverageAt(plane,strip,rig,therigb)); + // +} + +Float_t CaloFranzini::GetFullAverageAt(Int_t plane, Int_t strip, Float_t rig, Int_t therigb){ + // + Float_t minrig; + Float_t maxrig; + // + minrig = brigm->At(therigb); + maxrig = brigm->At(therigb+1); + // + Float_t minene = (*fqplmean[therigb])[plane][strip]; + Float_t maxene = (*fqplmean[therigb+1])[plane][strip]; + // + if ( maxrig == minrig ){ + printf("Unrecoverable ERROR! Matrix will be screwed up... \n"); + return(0.); + }; + Float_t b = log(maxene/minene)/(maxrig-minrig); + Float_t a = minene/exp(minrig*b); + Float_t ave = a*exp(b*rig); + if ( a == 0. || minene == 0. || ave != ave ){ + Float_t m = (maxene-minene)/(maxrig-minrig); + Float_t q = minene - m * minrig; + ave = rig * m + q; + }; + // + // ave += (44.-plane)*strip; + //if ( a == 0. ) ave = 0.; + if ( !(ave == ave) ) printf("a %f b %f ave %f maxene %f minene %f maxrig %f minrig %f \n",a,b,ave,maxene,minene,maxrig,minrig); // return(ave); // @@ -493,26 +1078,171 @@ // if ( rig < brigm->At(0) ){ if ( rig < brig->At(0) ){ - printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0)); + // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0)); }; therigb = 0; }; - if ( rig >= brigm->At(nbin-2) ){ + // if ( rig >= brigm->At(nbin-4) ){ // -2 + if ( rig >= brigm->At(nbin-2) ){ if ( rig >= brig->At(nbin-2) ) { - printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2)); + // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2)); }; + // therigb = nbin - 5;// -3 therigb = nbin - 3; }; // + if ( therigb < 2 ) therigb = 2; minrig = brigm->At(therigb); maxrig = brigm->At(therigb+1); + // printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex); // Float_t minene = (*hmat[therigb])[iindex][jindex]; Float_t maxene = (*hmat[therigb+1])[iindex][jindex]; + // printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave); + // +// Float_t a = 0.; +// Float_t b = 0.; +// Float_t ave = 0.; +// if ( minene == 0. ){ +// +// } else { +// b = log(maxene/minene)/(maxrig-minrig); +// a = minene/exp(minrig*b); +// ave = a*exp(b*rig); +// }; + // + Float_t m = (maxene-minene)/(maxrig-minrig); + Float_t q = minene - m * minrig; + Float_t ave = rig * m + q; + if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave); + // + // + return(ave); + // +} + +Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig){ + Int_t therigb = 0; + for (Int_t i = 0; i=brigm->At(i) && rig < brigm->At(i+1) ){ + therigb = i; + break; + }; + }; + // + if ( rig < brigm->At(0) ){ + if ( rig < brig->At(0) ){ + // printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brigm->At(0)); + }; + therigb = 0; + }; +// if ( rig >= brigm->At(nbin-4) ){ // -2 + if ( rig >= brigm->At(nbin-2) ){ + //if ( rig >= brigm->At(nbin-5) ){ + if ( rig >= brig->At(nbin-2) ) { + // printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brigm->At(nbin-2)); + }; + // therigb = nbin - 5;// -3 + // therigb = nbin - 5;// -3 + therigb = nbin - 3; + }; + // + if ( therigb < 2 ) therigb = 2; + if ( therigb > 13 ) therigb = 13; + // + return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb)); + // +} + + +Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb){ + // + return(this->GetFullHmatrixAt(iindex,jindex,rig,therigb,0)); + // +} + +Float_t CaloFranzini::GetFullHmatrixAt(Int_t iindex, Int_t jindex, Float_t rig, Int_t therigb, Int_t mtherig){ + // + // Int_t lofit = 12; + Int_t lofit = 3; + // + Float_t lowrig = 0.; + Float_t minrig; + Float_t maxrig; + // + if ( mtherig > lofit ) lowrig = brigm->At(therigb-1); + minrig = brigm->At(therigb); + maxrig = brigm->At(therigb+1); + //if ( therigb > 10 ) printf(" therigb %i minrig %f maxrig %f rig %f %i i %i j \n",therigb,minrig,maxrig,rig,iindex,jindex); + // + Float_t lowene = 0.; + if ( mtherig > lofit ) lowene = (*hfmat[therigb-1])[iindex][jindex]; + Float_t minene = (*hfmat[therigb])[iindex][jindex]; + Float_t maxene = (*hfmat[therigb+1])[iindex][jindex]; + // printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,a,b,rig,ave); + // + Float_t ave = 0.; + // + // Float_t a = 0.; + // Float_t b = 0.; + // Float_t ave = 0.; + // if ( minene == 0. ){ + // + // } else { + // b = log(maxene/minene)/(maxrig-minrig); + // a = minene/exp(minrig*b); + // ave = a*exp(b*rig); + // }; + // + if ( mtherig > lofit ){ + // + // FIT + // +// Float_t x[3]={lowrig,minrig,maxrig}; +// Float_t y[3]={lowene,minene,maxene}; +// // +// TGraph *tg= new TGraph(3,x,y); +// // +// // gStyle->SetLabelSize(0.04); +// // gStyle->SetNdivisions(510,"XY"); +// // TCanvas *c = new TCanvas(); +// // c->Draw(); +// // c->cd(); +// // tg->Draw("AP"); +// // +// TF1 *fun = new TF1("fun","pol2"); +// tg->Fit("fun","QNC"); +// // +// ave = fun->GetParameter(0) + rig * fun->GetParameter(1) + rig * rig * fun->GetParameter(2); +// // +// // printf(" therigb %i ave %f rig %f lowrig %f minrig %f maxrig %f lowene %f minene %f maxene %f \n",therigb,ave,rig,lowrig,minrig,maxrig,lowene,minene,maxene); +// // +// // +// // c->Modified(); +// // c->Update(); +// // gSystem->ProcessEvents(); +// // gSystem->Sleep(6000); +// // c->Close(); +// // gStyle->SetLabelSize(0); +// // gStyle->SetNdivisions(1,"XY"); +// // +// tg->Delete(); +// fun->Delete(); + Float_t mrigl = (lowrig+minrig)/2.; + Float_t mrigu = (minrig+maxrig)/2.; + Float_t menel = (lowene+minene)/2.; + Float_t meneu = (minene+maxene)/2.; + Float_t m = (meneu-menel)/(mrigu-mrigl); + Float_t q = menel - m * mrigl; + ave = rig * m + q; + // + } else { + Float_t m = (maxene-minene)/(maxrig-minrig); + Float_t q = minene - m * minrig; + ave = rig * m + q; + if ( debug ) printf(" therigb %i minrig %f maxrig %f minene %f maxene %f a %f b %f rig %f ave %f \n",therigb,minrig,maxrig,minene,maxene,m,q,rig,ave); + }; // - Float_t b = log(maxene/minene)/(maxrig-minrig); - Float_t a = minene/exp(minrig*b); - Float_t ave = a*exp(b*rig); // return(ave); // @@ -560,3 +1290,107 @@ gStyle->SetNdivisions(1,"XY"); // }; + +void CaloFranzini::DrawLongAverage(Int_t bin){ + // + TArrayF *ll = this->LoadLongAverage(brigm->At(bin)); + // + gStyle->SetLabelSize(0.04); + gStyle->SetNdivisions(510,"XY"); + // + TString hid = Form("cfralongvyvx"); + TCanvas *tcf = dynamic_cast(gDirectory->FindObject(hid)); + if ( tcf ){ + tcf->cd(); + } else { + tcf = new TCanvas(hid,hid); + }; + // + TString thid = Form("hfralongvyvx"); + TH1F *thf = dynamic_cast(gDirectory->FindObject(thid)); + if ( thf ) thf->Delete(); + thf = new TH1F(thid,thid,44,-0.5,43.5); + tcf->cd(); + Float_t qpl[43]; + for (Int_t st=0;st<43;st++){ + qpl[st] = ll->At(st); + printf("st %i qpl %f\n",st,qpl[st]); + }; + for (Int_t st=0;st<44;st++){ + if ( st == 37 ){ + thf->Fill(st,0.); + } else { + Int_t ss = st; + if ( st > 37 ) ss--; + thf->Fill(st,qpl[ss]); + }; + }; + thf->Draw(); + tcf->Modified(); + tcf->Update(); + // + gStyle->SetLabelSize(0); + gStyle->SetNdivisions(1,"XY"); + // +}; + +Int_t CaloFranzini::ConvertStrip(Int_t mstrip){ + // + Int_t lastrip = 0; +// // +// if ( mstrip < 50 ) lastrip = 0; +// // +// if ( mstrip >= 50 && mstrip < 64 ) lastrip = 1; +// // +// if ( mstrip >= 64 && mstrip < 70 ) lastrip = 2; +// // +// if ( mstrip >= 70 && mstrip < 75 ) lastrip = 3; +// // +// if ( mstrip >= 75 && mstrip < 84 ){ +// lastrip = (int)trunc((mstrip - 75)/3.) + 4; +// }; +// // +// if ( mstrip >= 84 && mstrip < 90 ){ +// lastrip = (int)trunc((mstrip - 84)/2.) + 7; +// }; +// // +// if ( mstrip >= 90 && mstrip < 101 ){ +// lastrip = mstrip - 90 + 10; +// }; +// // +// if ( mstrip >= 101 && mstrip < 107 ){ +// lastrip = (int)trunc((mstrip - 101)/2.) + 21; +// }; +// // +// // +// if ( mstrip >= 107 && mstrip < 116 ){ +// lastrip = (int)trunc((mstrip - 107)/3.) + 24; +// }; +// // +// if ( mstrip >= 116 && mstrip < 121 ) lastrip = 27; +// // +// if ( mstrip >= 121 && mstrip < 127 ) lastrip = 28; +// // +// if ( mstrip >= 127 && mstrip < 141 ) lastrip = 29; +// // +// if ( mstrip >= 141 ) lastrip = 30; +// // + // + if ( mstrip < 83 ) lastrip = 0; + // + if ( mstrip >= 83 && mstrip < 90 ) lastrip = 1; + // + if ( mstrip >= 90 && mstrip < 93 ) lastrip = 2; + // + if ( mstrip >= 93 && mstrip < 98 ){ + lastrip = mstrip - 93 + 3; + }; + // + if ( mstrip >= 98 && mstrip < 101 ) lastrip = 8; + // + if ( mstrip >= 101 && mstrip < 107 ) lastrip = 9; + // + if ( mstrip >= 107 ) lastrip = 10; + // + return(lastrip); +}