--- calo/flight/CaloProfile/src/CaloProfile.cpp 2007/11/09 09:11:24 1.1 +++ calo/flight/CaloProfile/src/CaloProfile.cpp 2008/09/22 20:08:25 1.2 @@ -1,4 +1,8 @@ #include +// +ClassImp(CaloLat); +ClassImp(CaloLong); +ClassImp(Calo2D); //////////////////////////////////////////////////////////////////////// /** @@ -60,7 +64,7 @@ /** * Default constructor */ -CaloLong::CaloLong(){ +Calo2D::Calo2D(){ Clear(); }; @@ -80,7 +84,7 @@ // }; -CaloLong::CaloLong(PamLevel2 *l2p){ +Calo2D::Calo2D(PamLevel2 *l2p){ // Clear(); // @@ -100,7 +104,7 @@ // }; -void CaloLong::Clear(){ +void Calo2D::Clear(){ // }; @@ -116,11 +120,11 @@ // }; -void CaloLong::Print(){ +void Calo2D::Print(){ // Process(); // - printf("==================== Calorimeter Lateral Profile =======================\n"); + printf("==================== Calorimeter 2D Profile =======================\n"); printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); // printf(" nx [number of X combination]:.. %i\n",nx); // printf(" ny [number of Y combination]:.. %i\n",ny); @@ -134,13 +138,12 @@ Draw(-1,-1); }; -void CaloLong::Draw(){ +void Calo2D::Draw(){ // Process(); Draw(-1); }; - void CaloLat::Draw(Int_t view,Int_t plane){ // Int_t minv = 0; @@ -203,16 +206,17 @@ // }; -void CaloLong::Draw(Int_t view){ +void Calo2D::Draw(Int_t plane){ // - Int_t minv = 0; - Int_t maxv = 0; + Int_t minp = 0; + Int_t maxp = 0; // - if ( view == -1 ){ - maxv = -1; + if ( plane == -1 ){ + minp = 0; + maxp = 23; } else { - minv = view; - maxv = view+1; + minp = plane; + maxp = plane+1; }; // Process(); @@ -220,59 +224,51 @@ gStyle->SetLabelSize(0.04); gStyle->SetNdivisions(510,"XY"); // - if ( maxv != -1 ){ - for (Int_t v=minv; v(gDirectory->FindObject(hid)); - if ( tc ){ -// tc->Clear(); - } else { - tc = new TCanvas(hid,hid); - }; - // - TString thid = Form("hlongv%i",v); - TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); - if ( th ) th->Delete(); -// th->Clear(); -// th->Reset(); -// } else { - th = new TH1F(thid,thid,22,-0.5,21.5); -// }; - tc->cd(); - // - for (Int_t st=0;st<22;st++){ - th->Fill(st,eplane[v][st]); - }; - th->Draw(); - tc->Modified(); - tc->Update(); - }; - } else { - // - TString hid = Form("clongvyvx"); - TCanvas *tc = dynamic_cast(gDirectory->FindObject(hid)); - if ( tc ){ - } else { - tc = new TCanvas(hid,hid); - }; - // - TString thid = Form("hlongvyvx"); - TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); - if ( th ) th->Delete(); - th = new TH1F(thid,thid,44,-0.5,43.5); - tc->cd(); - Int_t pp=0; - for (Int_t st=0;st<22;st++){ - for (Int_t v=1; v>=0;v--){ - // - th->Fill(pp,eplane[v][st]); - // - pp++; - }; - }; - th->Draw(); - tc->Modified(); - tc->Update(); + for (Int_t p=minp; p(gDirectory->FindObject(hid)); + if ( tc ){ + // tc->Clear(); + } else { + tc = new TCanvas(hid,hid); + }; + // + TString thid = Form("h2dp%i",p); + TH2F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + // th->Clear(); + // th->Reset(); + // } else { + Int_t minx = smax[p] - 10; + if ( minx < 0 ) minx = 0; + Int_t maxx = minx + 20; + if ( maxx > 95 ){ + maxx = 95; + minx = 75; + }; + Int_t miny = smay[p] - 10; + if ( miny < 0 ) miny = 0; + Int_t maxy = miny + 20; + if ( maxy > 95 ){ + maxy = 95; + miny = 75; + }; + th = new TH2F(thid,thid,20,(Float_t)minx-0.5,(Float_t)maxx-0.5,20,(Float_t)miny-0.5,(Float_t)maxy-0.5); + // th = new TH2F(thid,thid,96,-0.5,95.5,96,-0.5,95.5); + // }; + tc->cd(); + // + for (Int_t stx=minx;stxFill(stx,sty,estrip[p][stx][sty]); + }; + }; + gStyle->SetPalette(1); + // tc->SetLogz(); + // th->Draw("colbox"); + th->Draw("cont4"); + tc->Modified(); + tc->Update(); }; // gStyle->SetLabelSize(0); @@ -285,11 +281,12 @@ //delete this; }; -void CaloLong::Delete(){ +void Calo2D::Delete(){ Clear(); //delete this; }; + void CaloLat::Process(){ // if ( !L2 ){ @@ -335,7 +332,7 @@ // }; -void CaloLong::Process(){ +void Calo2D::Process(){ // if ( !L2 ){ printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); @@ -364,16 +361,565 @@ // // let's start // + Float_t es[2][22][96]; + memset(es,0, 4224*sizeof(Float_t)); + memset(estrip,0, 4224*sizeof(Float_t)); + Float_t mip1 = 0.; + Int_t view1 = 0; + Int_t plane1 = 0; + Int_t strip1 = 0; + // + for (Int_t i=0; iGetCaloLevel1()->istrip ; i++){ + mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1); + es[view1][plane1][strip1] = mip1; + }; + // + Int_t plane2 = 0; + Float_t emax[23]; + memset(emax,0,sizeof(Float_t)*23); + memset(smax,0,sizeof(Int_t)*23); + memset(smay,0,sizeof(Int_t)*23); + // + // + for (Int_t p=0; p < 23 ; p++){ + // + plane1 = p-1; + plane2 = p; + // + if ( p == 0 ){ + plane1 = -1; + plane2 = 0; + }; + if ( p == 22 ){ + plane1 = 21; + plane2 = -1; + }; + // + for (Int_t s=0; s < 96 ; s++){ // x + for (Int_t ss=0; ss < 96 ; ss++){ // y + if ( p == 0 ){ + estrip[p][s][ss] += es[1][plane2][ss]; + if ( (es[1][plane2][ss]) > emax[p] ){ + smax[p] = 45; + smay[p] = ss; + emax[p] = es[1][plane2][ss] ; + }; + }; + if ( p > 0 && p < 22 ){ + estrip[p][s][ss] += es[0][plane1][s] + es[1][plane2][ss]; + if ( (es[0][plane1][s] + es[1][plane2][ss]) > emax[p] ){ + smax[p] = s; + smay[p] = ss; + emax[p] = es[0][plane1][s] + es[1][plane2][ss] ; + }; + }; + if ( p == 22 ){ + estrip[p][s][ss] += es[0][plane1][s]; + if ( (es[1][plane2][s]) > emax[p] ){ + smax[p] = s; + smay[p] = 45; + emax[p] = es[1][plane2][s] ; + }; + }; + }; + }; + // + }; + // + if ( debug ) this->Print(); + if ( debug ) printf(" exit \n"); + // +}; + + +/** + * Default constructor + */ +CaloLong::CaloLong(){ + Clear(); +}; + +CaloLong::CaloLong(PamLevel2 *l2p){ + // + Clear(); + // + L2 = l2p; + // + if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n"); + // + OBT = 0; + PKT = 0; + atime = 0; + // + sel = true; + cont = false; + N = 0; + NC = 22; + mask18b = -1; + // + no18x = true; + debug = false; + // +}; + +void CaloLong::Clear(){ + // memset(eplane,0, 2*22*sizeof(Float_t)); // + chi2 = 0.; + ndf = 0.; + E0 = 0.; + a = 0.; + b = 0.; + errE0 = 0.; + erra = 0.; + errb = 0.; + etmax = 0.; + asymm = 0.; + fitresult = 0; + // + X0pl = 0.76; + // +}; + +void CaloLong::Print(){ + // + Process(); + // + printf("==================== Calorimeter Longitudinal Profile =======================\n"); + printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime); + printf(" fitresult:.. %i\n",fitresult); + printf(" chi2 :.. %f\n",chi2); + printf(" ndf :.. %f\n",ndf); + printf(" nchi2 :.. %f\n",chi2/ndf); + printf(" E0 :.. %f\n",E0); + printf(" E0/260. :.. %f\n",E0/260.); + printf(" a :.. %f\n",a); + printf(" b :.. %f\n",b); + printf(" errE0 :.. %f\n",errE0); + printf(" erra :.. %f\n",erra); + printf(" errb :.. %f\n",errb); + printf(" asymm :.. %f\n",asymm); + printf(" tmax :.. %f\n",((a-1.)/b)); + printf(" etmax :.. %f\n",etmax); + printf(" X0pl :.. %f\n",X0pl); + printf("========================================================================\n"); + // +}; + +void CaloLong::SetNoWpreSampler(Int_t n){ + // + if ( NC+n <= 22 && NC+n >= 0 ){ + 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 = 22; + N = 0; + }; +} + +void CaloLong::SetNoWcalo(Int_t n){ + if ( N+n <= 22 && N+n >= 0 ){ + NC = n; + } 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 = 22; + N = 0; + }; +} + +void CaloLong::Process(){ + // + if ( !L2 ){ + printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); + printf(" ERROR: CaloHough variables not filled \n"); + return; + }; + // + Bool_t newentry = false; + // + if ( L2->IsORB() ){ + if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){ + newentry = true; + OBT = L2->GetOrbitalInfo()->OBT; + PKT = L2->GetOrbitalInfo()->pkt_num; + atime = L2->GetOrbitalInfo()->absTime; + }; + } else { + newentry = true; + }; + // + if ( !newentry ) return; + // + if ( debug ) printf(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime); + // + Clear(); + // + // let's start + // + if ( cont ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18+N) ){ + mask18b = 18 + N; + break; + }; + }; + }; + // + if ( sel ){ + for (Int_t i=0; i<22; i++){ + if ( i == (18-N) ){ + mask18b = 18 - N; + break; + }; + }; + }; + // + // if ( mask18b == 18 ) mask18b = -1; + // + Int_t view = 0; + Int_t plane = 0; + Int_t strip = 0; + Float_t mip = 0.; + for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){ + mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip); + eplane[view][plane] += mip; + }; + // + // inclination factor (steal from Daniele's code) + // + Float_t ytgx = 0; + Float_t ytgy = 0; + ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0]; + ytgy = 0.76 * L2->GetCaloLevel2()->tany[0]; + X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) ); + // + // Find experimental plane of maximum + // + Int_t pmax = 0; + Int_t vmax = 0; + Float_t emax = 0.; for (Int_t v=0; v<2; v++){ for (Int_t i=0; i<22; i++){ - eplane[v][i] = L2->GetCaloLevel1()->qtotpl(v,i); + if ( eplane[v][i] > emax ){ + emax = eplane[v][i]; + vmax = v; + pmax = i; + }; }; }; // + // + // + if ( vmax == 0 ) pmax++; + etmax = pmax * X0pl; + // if ( debug ) this->Print(); if ( debug ) printf(" exit \n"); // }; + +Double_t ccurve(Double_t *ti,Double_t *par){ + // + Double_t t = *ti; + Double_t cE0 = par[0]; + Double_t ca = par[1]; + Double_t cb = par[2]; + Double_t gammaa = TMath::Gamma(ca); + // + Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa; + // + return value; + // +} + +void CaloLong::Fit(){ + this->Fit(false); +}; + +void CaloLong::Fit(Bool_t draw){ + // + Process(); + // + // + if ( !L2 ){ + printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n"); + printf(" ERROR: CaloHough variables not filled \n"); + return; + }; + // + Bool_t newentry = false; + // + if ( L2->IsORB() ){ + if ( L2->GetOrbitalInfo()->pkt_num != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){ + newentry = true; + fOBT = L2->GetOrbitalInfo()->OBT; + fPKT = L2->GetOrbitalInfo()->pkt_num; + fatime = L2->GetOrbitalInfo()->absTime; + }; + } else { + newentry = true; + }; + // + if ( !newentry ) return; + // + if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime); + // + if ( draw ){ + gStyle->SetLabelSize(0.04); + gStyle->SetNdivisions(510,"XY"); + }; + // + TString hid = Form("clongfit"); + TCanvas *tc = dynamic_cast(gDirectory->FindObject(hid)); + // if ( tc ) tc->Delete(); + // if ( tc ) tc->Close(); + if ( !tc && draw ){ + tc = new TCanvas(hid,hid); + } else { + if ( tc ) tc->cd(); + }; + // + TString thid = Form("hlongfit"); + TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + Float_t xpos = 0.; + Float_t enemip = 0.; + Float_t xmax = NC * X0pl + 0.2; + // th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax); + th = new TH1F(thid,thid,100,-0.2,xmax); + // + for (Int_t st=N;st<(N+NC);st++){ + enemip = 0.; + xpos = (st - N) * X0pl; + if ( st > N && st < (N+NC-1) ){ + if ( no18x && ( st == 18+1 || st == mask18b+1 )){ + enemip = 2. * eplane[1][st]; + } else { + enemip = eplane[0][st-1] + eplane[1][st]; + }; + } else { + if ( st == N ) enemip = 2. * eplane[1][st]; + if ( st == (N+NC-1) ) enemip = 2. * eplane[0][st]; + }; + // + if ( enemip > 0. ){ + th->Fill(xpos,enemip); + if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip); + }; + // + // for (Int_t v=1; v>=0;v--)// { +// // +// if ( v == 1 ){ +// xpos = (st - N) * X0pl; +// } else { +// xpos = (st + 1 - N) * X0pl; +// }; +// // +// if ( no18x && st == 18 && v == 0 ){ +// // skip plane 18x +// } else { +// if ( v == 1 && st == mask18b ){ +// // emulate plane 18x +// } else { +// if ( eplane[v][st] > 0. ){ +// th->Fill(xpos,eplane[v][st]); +// if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]); +// }; +// }; +// }; +// // +// }; + }; + // + TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3); + E0 = L2->GetCaloLevel2()->qtot; + a = 5.; + b = 0.5; + if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b); + lfit->SetParameters(E0,a,b); + // lfit->SetParLimits(0,0.,1000.); + // lfit->SetParLimits(1,-1.,80.); + // lfit->SetParLimits(2,-1.,10.); + TString optio; + if ( debug ){ + // optio = "MERBOV"; + // optio = "MEROV"; + // optio = "EROV"; + optio = "RNOV"; + if ( draw ) optio = "ROV"; + } else { + // optio = "MERNOQ"; + // optio = "ERNOQ"; + optio = "RNOQ"; + if ( draw ) optio = "ROQ"; + }; + // + if ( debug ) printf(" OK, start the fitting procedure...\n"); + // + fitresult = th->Fit("lfit",optio); + // + if ( debug ) printf(" the fit is done! result: %i \n",fitresult); + // + E0 = lfit->GetParameter(0); + a = lfit->GetParameter(1); + b = lfit->GetParameter(2); + errE0 = lfit->GetParError(0); + erra = lfit->GetParError(1); + errb = lfit->GetParError(2); + chi2 = lfit->GetChisquare(); + ndf = lfit->GetNDF(); + Float_t tmax = 0.; + if ( debug ) printf(" Parameters are retrieved \n"); + if ( b != 0 ) tmax = (a - 1.)/b; + // + if ( fitresult != 0 ){ + if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n"); + asymm = -1.; + } else { + Int_t npp = 1000; + double *xp=new double[npp]; + double *wp=new double[npp]; + lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12); + Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax); + // Float_t imax = lfit->Integral(0.,tmax); + if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax); + Int_t np = 1000; + double *x=new double[np]; + double *w=new double[np]; + lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12); + Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax); + delete x; + delete w; + delete xp; + delete wp; + // Float_t i10max = lfit->Integral(0.,10.*tmax); + if ( debug ) printf(" Integral: %f \n",i10max); + // + if ( i10max != imax ){ + asymm = imax / (i10max-imax); + } else { + if ( debug ) printf(" i10max == imax, asymm undefined\n"); + asymm = -2.; + }; + //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax)); + if ( debug ) printf(" Asymmetry has been calculated \n"); + }; + // + if ( draw ){ + // + tc->cd(); + // gStyle->SetOptStat(11111); + tc->SetTitle(); + th->SetTitle(""); + th->SetName(""); + th->SetMarkerStyle(20); + // axis titles + th->SetXTitle("Depth [X0]"); + th->SetYTitle("Energy [MIP]"); + th->DrawCopy("Perror"); + lfit->Draw("same"); + tc->Modified(); + tc->Update(); + // + gStyle->SetLabelSize(0); + gStyle->SetNdivisions(1,"XY"); + // + }; + // + delete lfit; + // +}; + +void CaloLong::Draw(){ + // + Process(); + Draw(-1); +}; + +void CaloLong::Draw(Int_t view){ + // + Int_t minv = 0; + Int_t maxv = 0; + // + if ( view == -1 ){ + maxv = -1; + } else { + minv = view; + maxv = view+1; + }; + // + Process(); + // + gStyle->SetLabelSize(0.04); + gStyle->SetNdivisions(510,"XY"); + // + if ( maxv != -1 ){ + for (Int_t v=minv; v(gDirectory->FindObject(hid)); + if ( tc ){ + // tc->Clear(); + } else { + tc = new TCanvas(hid,hid); + }; + // + TString thid = Form("hlongv%i",v); + TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + // th->Clear(); + // th->Reset(); + // } else { + th = new TH1F(thid,thid,22,-0.5,21.5); + // }; + tc->cd(); + // + for (Int_t st=0;st<22;st++){ + th->Fill(st,eplane[v][st]); + }; + th->Draw(); + tc->Modified(); + tc->Update(); + }; + } else { + // + TString hid = Form("clongvyvx"); + TCanvas *tc = dynamic_cast(gDirectory->FindObject(hid)); + if ( tc ){ + } else { + tc = new TCanvas(hid,hid); + }; + // + TString thid = Form("hlongvyvx"); + TH1F *th = dynamic_cast(gDirectory->FindObject(thid)); + if ( th ) th->Delete(); + th = new TH1F(thid,thid,44,-0.5,43.5); + tc->cd(); + Int_t pp=0; + for (Int_t st=0;st<22;st++){ + for (Int_t v=1; v>=0;v--){ + // + th->Fill(pp,eplane[v][st]); + // + pp++; + }; + }; + th->Draw(); + tc->Modified(); + tc->Update(); + }; + // + gStyle->SetLabelSize(0); + gStyle->SetNdivisions(1,"XY"); + // +}; + +void CaloLong::Delete(){ + Clear(); + //delete this; +};