--- calo/flight/CaloFranzini/src/CaloFranzini.cpp 2007/12/04 12:42:54 1.2 +++ calo/flight/CaloFranzini/src/CaloFranzini.cpp 2007/12/18 09:55:07 1.4 @@ -18,6 +18,7 @@ // file = NULL; brig = NULL; + brigm = NULL; nbin = 0; // L2 = l2p; @@ -31,6 +32,13 @@ dofull = false; sntr = 0; // + sel = true; + cont = false; + N = 0; + NC = 43; + // + mask18b = -1; + // Clear(); // } @@ -44,7 +52,7 @@ fulltzeta = -100.; degfre = 0; memset(estrip, 0, 2*22*96*sizeof(Float_t)); - memset(qplane, 0, 2*22*sizeof(Float_t)); + memset(qplane, 0, 43*sizeof(Float_t)); // } @@ -55,8 +63,11 @@ 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(" longtzeta :.. %f\n",longtzeta); - printf(" fulltzeta :.. %f\n",fulltzeta); + printf(" longtzeta normalized :.. %f\n",this->GetNormLongTZeta()); + // printf(" fulltzeta :.. %f\n",fulltzeta); + // printf(" longtzeta normalized :.. %f\n",this->GetNormFullTZeta()); printf("========================================================================\n"); // } @@ -69,6 +80,30 @@ // } +void CaloFranzini::SetNoWpreSampler(Int_t n){ + if ( NC+n < 23 ){ + N = n; + } else { + printf(" ERROR! Calorimeter is made of 22 W planes\n"); + printf(" you are giving N presampler = %i and N calo = %i \n",n,NC); + printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n"); + NC = 43; + N = 0; + }; +} + +void CaloFranzini::SetNoWcalo(Int_t n){ + if ( N+n < 23 ){ + NC = n*2; + if ( NC >37 ) NC--; + } else { + printf(" ERROR! Calorimeter is made of 22 W planes\n"); + printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n); + printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n"); + NC = 43; + N = 0; + }; +} void CaloFranzini::Process(){ this->Process(0); @@ -124,8 +159,12 @@ mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); // estrip[view][plane][strip] = mip; - nplane = 1 - view + 2 * plane; - qplane[nplane] += 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; // }; // @@ -133,25 +172,51 @@ // if ( dolong ){ // - degfre = 44; + 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; + // + Int_t dgf = 43; // for (Int_t i=0; i < 22; i++){ if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][1] < 0 ){ - degfre = 2 * i; + dgf = 2 * i; break; }; if ( L2->GetTrack(itr)->GetCaloTrack()->tibar[i][0] < 0 ){ - degfre = 1 + 2 * i; + dgf = 1 + 2 * i; break; }; }; // + if ( dgf < 43 && dgf > 37 ) dgf--; + // + degfre = TMath::Min(dgf,NC); + // if ( degfre > 0 ){ - TArrayF *qplmean = this->LoadLongAverage(rig); - TMatrixD *lmatrix = this->LoadCovarianceMatrix(rig); for (Int_t i = 0; i < degfre; i++){ - for (Int_t j = 0; j < degfre; j++){ - longtzeta += (qplane[i]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j)); + if ( i != mask18b ){ + for (Int_t j = 0; j < degfre; j++){ + if ( j != mask18b ){ + longtzeta += (qplane[i]-this->GetAverageAt(i,rig)) * this->GetHmatrixAt(i,j,rig) * (qplane[j]-this->GetAverageAt(j,rig)); + }; + }; }; }; }; @@ -186,6 +251,7 @@ // if ( !file || file->IsZombie() ){ file = new TFile(matrixfile.Data(),"RECREATE"); + printf(" Create file %s \n",file->GetName()); } else { printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data()); return(false); @@ -195,26 +261,53 @@ // } +Bool_t CaloFranzini::UpdateMatrixFile(TString matrixfile){ + // + file = new TFile(matrixfile.Data(),"UPDATE"); + // + if ( !file || file->IsZombie() ){ + printf(" ERROR: file %s already existing!\n Choose another name or delete the old file\n",matrixfile.Data()); + return(false); + }; + // + return(true); + // +} + void CaloFranzini::WriteNumBin(Int_t numbin){ + file->cd(); TArrayI nbi(1, &numbin); file->WriteObject(&nbi, "nbinenergy"); } void CaloFranzini::WriteRigBin(TArrayF *rigbin){ - file->WriteObject(&rigbin, "binrig"); + file->cd(); + // rigbin->Write("binrig"); + file->WriteObject(&(*rigbin), "binrig"); } void CaloFranzini::WriteLongMean(TArrayF *qpl, Int_t bin){ + file->cd(); TString name = Form("qplmeann%i",bin); - file->WriteObject(&qpl,name.Data()); + file->WriteObject(&(*qpl),name.Data()); } -void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){ +void CaloFranzini::WriteInvertedLongMatrix(TMatrixD mat, Int_t bin){ + file->cd(); TString name = Form("matrixn%i",bin); + // mat.Write(name.Data()); + file->WriteObject(&mat,name.Data()); +} + +void CaloFranzini::WriteLongMatrix(TMatrixD *mat, Int_t bin){ + file->cd(); + TString name = Form("origmatrixn%i",bin); + // mat.Write(name.Data()); file->WriteObject(&(*mat),name.Data()); } void CaloFranzini::CloseMatrixFile(){ + file->cd(); file->Close(); } @@ -235,76 +328,235 @@ return(false); }; // + if ( !this->LoadBin() ){ + printf(" %s \n",matrixfile.Data()); + return(false); + }; + if ( !this->LoadMatrices() ){ + printf(" %s \n",matrixfile.Data()); + return(false); + }; + // + // + return(true); + // +} + +Bool_t CaloFranzini::LoadBin(){ + // TArrayI *numbin = (TArrayI*)file->Get("nbinenergy"); if ( !numbin ){ - printf(" ERROR: cannot read number of bins from file %s \n",matrixfile.Data()); + 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 %s) \n",matrixfile.Data()); + printf(" ERROR: cannot work with 0 energy bins from file "); return(false); }; // brig = (TArrayF*)file->Get("binrig"); if ( !brig ){ - printf(" ERROR: cannot read rigidity binning from file %s \n",matrixfile.Data()); + printf(" ERROR: cannot read rigidity binning from file "); + return(false); + }; + // + 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*)file->Get(name.Data()); + if ( !qplmean[i] ){ + printf(" ERROR: cannot read average from file "); + return(false); + }; + }; + // return(true); +} + +Bool_t CaloFranzini::LoadMatrices(){ + // + for (Int_t i=0;i<17;i++){ + TString name1 = Form("matrixn%i",i); + hmat[i] = (TMatrixD*)file->Get(name1.Data()); + }; // + return(true); } TMatrixD *CaloFranzini::LoadCovarianceMatrix(Float_t rig){ // - TString name; + Int_t mv = 0; for (Int_t i = 0; i=brig->At(i) && rig < brig->At(i+1) ){ - name = Form("matrixn%i",i); + mv = i; break; }; }; if ( rig < brig->At(0) ){ printf(" WARNING: Event with rigidity lower than the first covariance matrix bin (rig = %f, lower limit = %f)\n",rig,brig->At(0)); - name = "matrixn0"; - printf(" Using matrix %s \n",name.Data()); + mv = 0; }; - if ( rig >= brig->At(nbin) ){ - printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin)); - name = Form("matrixn%i",nbin-1); - printf(" Using matrix %s \n",name.Data()); + if ( rig >= brig->At(nbin-1) ){ + printf(" WARNING: Event with rigidity higher than the last covariance matrix bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1)); + mv = nbin-2; }; // - TMatrixD *matrix = (TMatrixD*)file->Get(name.Data()); - // - return(matrix); + return(hmat[mv]); // } TArrayF *CaloFranzini::LoadLongAverage(Float_t rig){ // - TString name; + Int_t mv=0; for (Int_t i = 0; i=brig->At(i) && rig < brig->At(i+1) ){ - name = Form("qplmeann%i",i); + mv = i; break; }; }; if ( rig < brig->At(0) ){ printf(" WARNING: Event with rigidity lower than the first qplmean bin (rig = %f, lower limit = %f)\n",rig,brig->At(0)); - name = "qplmeann0"; - printf(" Using qplmean %s \n",name.Data()); + mv = 0; + }; + if ( rig >= brig->At(nbin-1) ){ + printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin-1)); + mv=nbin-2; + }; + // + return(qplmean[mv]); + // +} + +Float_t CaloFranzini::GetAverageAt(Int_t plane, Float_t rig){ + // + Int_t therigb = 0; + for (Int_t i = 0; i=brigm->At(i) && rig < brigm->At(i+1) ){ + therigb = i; + break; + }; + }; + // + Float_t minrig; + Float_t maxrig; + // + // + 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; + }; + // + minrig = brigm->At(therigb); + maxrig = brigm->At(therigb+1); + // + Float_t minene = (*qplmean[therigb])[plane]; + Float_t maxene = (*qplmean[therigb+1])[plane]; + // + if ( maxrig == minrig ){ + printf("Unrecoverable ERROR! Matrix will be screwed up... \n"); + return(0.); }; - if ( rig >= brig->At(nbin) ){ - printf(" WARNING: Event with rigidity higher than the last qplmean bin (rig = %f, upper limit = %f)\n",rig,brig->At(nbin)); - name = Form("qplmeann%i",nbin-1); - printf(" Using qplmean %s \n",name.Data()); + Float_t b = log(maxene/minene)/(maxrig-minrig); + Float_t a = minene/exp(minrig*b); + Float_t ave = a*exp(b*rig); + // + return(ave); + // +} + +Float_t CaloFranzini::GetHmatrixAt(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; + }; }; // - TArrayF *qplmean = (TArrayF*)file->Get(name.Data()); + Float_t minrig; + Float_t maxrig; + // + 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; + }; + // + minrig = brigm->At(therigb); + maxrig = brigm->At(therigb+1); + // + Float_t minene = (*hmat[therigb])[iindex][jindex]; + Float_t maxene = (*hmat[therigb+1])[iindex][jindex]; // - return(qplmean); + Float_t b = log(maxene/minene)/(maxrig-minrig); + Float_t a = minene/exp(minrig*b); + Float_t ave = a*exp(b*rig); + // + return(ave); // } + +void CaloFranzini::DrawLongAverage(Float_t rig){ + // + TArrayF *ll = this->LoadLongAverage(rig); + // + 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"); + // +};