--- calo/flight/CaloFranzini/src/CaloFranzini.cpp 2007/12/04 12:42:54 1.2 +++ calo/flight/CaloFranzini/src/CaloFranzini.cpp 2007/12/13 17:08:11 1.3 @@ -31,6 +31,13 @@ dofull = false; sntr = 0; // + sel = true; + cont = false; + N = 0; + NC = 43; + // + mask18b = -1; + // Clear(); // } @@ -44,7 +51,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 +62,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 +79,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 +158,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 +171,53 @@ // 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]-qplmean->At(i)) * (*lmatrix)[i][j] * (qplane[j]-qplmean->At(j)); + }; + }; }; }; }; @@ -186,6 +252,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); @@ -196,25 +263,39 @@ } 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(); } @@ -270,9 +351,9 @@ name = "matrixn0"; printf(" Using matrix %s \n",name.Data()); }; - 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); + 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)); + name = Form("matrixn%i",nbin-2); printf(" Using matrix %s \n",name.Data()); }; // @@ -297,9 +378,9 @@ name = "qplmeann0"; printf(" Using qplmean %s \n",name.Data()); }; - 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); + 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)); + name = Form("qplmeann%i",nbin-2); printf(" Using qplmean %s \n",name.Data()); }; // @@ -308,3 +389,48 @@ return(qplmean); // } + + + +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(); + // Int_t pp=0; + Float_t qpl[43]; + for (Int_t st=0;st<43;st++){ + qpl[st] = ll->At(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"); + // +};