--- quicklook/tof/src/TofScan.cpp 2006/05/11 10:07:39 1.1.1.1 +++ quicklook/tof/src/TofScan.cpp 2007/10/08 12:28:54 1.4 @@ -4,6 +4,11 @@ * Version 1.2 * Modified by G.De Rosa * Date 27 Apr 2006 +* Modified by G.De Rosa +* Date 03 Jul 2006 +* Modified by W. Menn to select helium particles for PMT gain check +* Date 09 Aug 2007 +* Last version 08 Oct 2007 * * Description: * Describe the performance of the TOF. @@ -14,6 +19,10 @@ * TString format - the format which will be used fo rsave the produced images (Default = "gif") */ +#include +#include +#include +#include #include #include #include @@ -25,8 +34,9 @@ #include #include - #include +#include + using namespace std; void TofScan(TString base, TString outDirectory = "", TString format = ""){ @@ -52,10 +62,6 @@ file->Close(); return; } - - TH2S *tdcHisto = new TH2S("tdcHisto", sst.str().c_str(), 409, 0, 4096, 48, 0, 48); - TH2S *adcHisto = new TH2S("adcHisto", sst.str().c_str(), 409, 0, 4096, 48, 0, 48); - /* * Array to convert hdc/adc to the real Photomultiplier @@ -94,69 +100,361 @@ "S32_1A", "S32_1B", "S32_2A", "S32_2B", "S32_3A", "S32_3B" }; - int i = 0; + const Int_t nh = 48; + TH1F *htdc[nh]; + TH1F *hadc[nh]; + + TObjArray *hhtdc = new TObjArray(nh); + TObjArray *hhadc = new TObjArray(nh); + char tdcname[48]=""; + char adcname[48]=""; + + char htitle[50]; + TH1F *adche[48]; + for(int i=0;i<48;i++) { + sprintf(htitle, "adche_%d",(i+1)); + adche[i] = new TH1F(htitle,htitle,100,0.,1500.); + } + + + Float_t adca[48]; // vector with adc values according to "ind"=pmt_id + Float_t tdca[48]; // the same for tdc + int j = 0; int k = 0; int z = 0; int ch = 0; int hb = 0; - - while(i < nevents){ + int ind =0; + + int heevent =0; + +// upper and lower limits for the helium selection +Float_t A_l[24]={200,190,300,210,220,200,210,60, 60, 120,220,120,160,50, 300,200,120,250,350,300,350,250,280,300}; +Float_t A_h[24]={550,490,800,600,650,600,600,260,200,380,620,380,550,200,850,560,400,750,900,800,880,800,750,800}; + +// The k1 constants for the beta calculation, only for S1-S3 +// k2 constant is taken to be the standard 2D/c +Float_t k1[72] = {50,59.3296,28.4328,-26.0818,5.91253,-19.588,-9.26316,24.7544,2.32465, +-50.5058,-15.3195,-39.1443,-91.2546,-58.6243,-84.5641,-63.1516,-32.2091, +-58.3358,13.8084,45.5322,33.2416,-11.5313,51.3271,75,-14.1141, +42.8466,15.1794,-63.6672,-6.07739,-32.164,-41.771,10.5274,-9.46096, +-81.7404,-28.783,-52.7167,-127.394,-69.6166,-93.4655,-98.9543,-42.863, +-67.8244,-19.3238,31.1221,8.7319,-43.1627,5.55573,-14.4078,-83.4466, +-47.4647,-77.8379,-108.222,-75.986,-101.297,-96.0205,-63.1881,-90.1372, +-22.7347,8.31409,-19.6912,-7.49008,23.6979,-1.66677,1.81556,34.4668, +6.23693,-100,-59.5861,-90.9159,-141.639,-89.2521,-112.881} ; + +//------------------------------------------------------------------- + + + for (int i=0; i < nevents; i++){ + TofBr->GetEntry(i); + k = 0; while (k < 24){ j = 0; while (j < 2){ ch = tof[2*j][k] - 1; hb = tof[2*j + 1][k] - 1; - tdcHisto->Fill(tofEvent->tdc[ch][hb], 2*k + j); - adcHisto->Fill(tofEvent->adc[ch][hb], 2*k + j); + ind = 2*k + j; + + if(i==0){ + sprintf(tdcname,"TDChist%4.4d",ind); + sprintf(adcname,"ADChist%4.4d",ind); + + htdc[ind] = new TH1F(tdcname,tdcname,409,0,4096); + hadc[ind] = new TH1F(adcname,adcname,409,0,4096); + + hhtdc->Add(htdc[ind]); + hhadc->Add(hadc[ind]); + } + + htdc[ind]->Fill(tofEvent->tdc[ch][hb]); + hadc[ind]->Fill(tofEvent->adc[ch][hb]); + tdca[ind]=tofEvent->tdc[ch][hb]; + adca[ind]=tofEvent->adc[ch][hb]; j++; } k++; } - i++; - } - - TH1D *tpADC = 0; - TH1D *tpTDC = 0; - TH1D *tpADCall = 0; - TH1D *tpTDCall = 0; - + +//============ calculate beta and select helium ==================== + +// find hitted paddle by looking for ADC values on both sides +// since we looking for helium this gives decent results + +Int_t tof11_i,tof12_i,tof21_i,tof22_i,tof31_i,tof32_i; +Float_t a1,a2; +Int_t jj; + +// reset values +tof11_i = -1; +tof12_i = -1; +tof21_i = -1; +tof22_i = -1; +tof31_i = -1; +tof32_i = -1; + + for(jj=0; jj<8; jj++){ + a1 = adca[2*jj]; + a2 = adca[2*jj+1]; + if ((a1 < 3000) && (a2 < 3000)) tof11_i = jj; + } + for(jj=0; jj<6; jj++){ + a1 = adca[16+2*jj]; + a2 = adca[16+2*jj+1]; + if ((a1 < 3000) && (a2 < 3000)) tof12_i = jj; + } + for(jj=0; jj<2; jj++){ + a1 = adca[28+2*jj]; + a2 = adca[28+2*jj+1]; + if ((a1 < 3000) && (a2 < 3000)) tof21_i = jj; + } + for(jj=0; jj<2; jj++){ + a1 = adca[32+2*jj]; + a2 = adca[32+2*jj+1]; + if ((a1 < 3000) && (a2 < 3000)) tof22_i = jj; + } + for(jj=0; jj<3; jj++){ + a1 = adca[36+2*jj]; + a2 = adca[36+2*jj+1]; + if ((a1 < 3000) && (a2 < 3000)) tof31_i = jj; + } + for(jj=0; jj<3; jj++){ + a1 = adca[42+2*jj]; + a2 = adca[42+2*jj+1]; + if ((a1 < 3000) && (a2 < 3000)) tof32_i = jj; + } + + +//---------------------------------------------------------------- + +Float_t zin[6] = {53.74, 53.04, 23.94, 23.44, -23.49, -24.34}; +Float_t c1,c2,xhelp,xhelp1,xhelp2,ds,dist,F; +Float_t sw,sxw,beta_mean_tof,w_i; +Float_t theta,x1,x2,y1,y2,dx,dy,dr; +Int_t ihelp; +Int_t ipmt[4]; +Float_t time[4]; +Float_t beta1[4]; + +// Only use events with: S11 and S12 and S31 and S32 + +if ( (tof11_i>-1) && (tof12_i>-1) && (tof31_i>-1) && (tof32_i>-1) ) { + +// calculate zenith angle theta using the locations of the hitted paddles + + + Float_t tof11_x[8] = {-17.85,-12.75,-7.65,-2.55,2.55,7.65,12.75,17.85}; + Float_t tof12_y[6] = { -13.75,-8.25,-2.75,2.75,8.25,13.75}; +// Float_t tof21_y[2] = { 3.75,-3.75}; +// Float_t tof22_x[2] = { -4.5,4.5}; + Float_t tof31_x[3] = { -6.0,0.,6.0}; + Float_t tof32_y[3] = { -5.0,0.0,5.0}; + + // S11 8 paddles 33.0 x 5.1 cm + // S12 6 paddles 40.8 x 5.5 cm + // S21 2 paddles 18.0 x 7.5 cm + // S22 2 paddles 15.0 x 9.0 cm + // S31 3 paddles 15.0 x 6.0 cm + // S32 3 paddles 18.0 x 5.0 cm + + x1 = 0.; + x2 = 0.; + y1 = 0.; + y2 = 0.; + + x1 = tof11_x[tof11_i] ; + y1 = tof12_y[tof12_i] ; + x2 = tof31_x[tof31_i] ; + y2 = tof32_y[tof32_i] ; + + theta=0.; + dx=0.; + dy=0.; + dr=0.; + + dx = x1-x2; + dy = y1-y2; + dr = sqrt(dx*dx+dy*dy); + theta = atan(dr/77.5); + + + beta_mean_tof=100.; + + for (Int_t jj=0; jj< 4; jj++) beta1[jj] = 100. ; + + +//---------------------------------------------------------------- +//--------- S1 - S3 --------------------------------------------- +//---------------------------------------------------------------- + +//--------- S11 - S31 ------------------------------------------- + +if ((tof11_i>-1)&&(tof31_i>-1)) { + +dist = zin[0] - zin[4]; +c2 = (2.*0.01*dist)/(3.E08*50.E-12); +F = 1./cos(theta); + + ipmt[0] = (tof11_i)*2; + ipmt[1] = (tof11_i)*2+1; + ipmt[2] = 36+(tof31_i)*2; + ipmt[3] = 36+(tof31_i)*2+1; + + for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ; + + if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) { + xhelp1 = time[0] + time[1] ; + xhelp2 = time[2] + time[3] ; + ds = xhelp1-xhelp2; + ihelp=0+(tof11_i)*3+tof31_i ; + c1 = k1[ihelp] ; + beta1[0] = c2*F/(ds-c1); +} +} + +//--------- S11 - S32 ------------------------------------------- + +if ((tof11_i>-1)&&(tof32_i>-1)) { + +dist = zin[0] - zin[5]; +F = 1./cos(theta); +c2 = (2.*0.01*dist)/(3.E08*50.E-12); + + ipmt[0] = (tof11_i)*2; + ipmt[1] = (tof11_i)*2+1; + ipmt[2] = 42+(tof32_i)*2; + ipmt[3] = 42+(tof32_i)*2+1; + + for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ; + + if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) { + xhelp1 = time[0] + time[1] ; + xhelp2 = time[2] + time[3] ; + ds = xhelp1-xhelp2; + ihelp=24+(tof11_i)*3+tof32_i ; + c1 = k1[ihelp] ; + beta1[1] = c2*F/(ds-c1); +} +} + +//--------- S12 - S31 ------------------------------------------- + +if ((tof12_i>-1)&&(tof31_i>-1)) { + +dist = zin[1] - zin[4]; +F = 1./cos(theta); +c2 = (2.*0.01*dist)/(3.E08*50.E-12); + + ipmt[0] = 16+(tof12_i)*2; + ipmt[1] = 16+(tof12_i)*2+1; + ipmt[2] = 36+(tof31_i)*2; + ipmt[3] = 36+(tof31_i)*2+1; + + for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ; + + if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) { + xhelp1 = time[0] + time[1] ; + xhelp2 = time[2] + time[3] ; + ds = xhelp1-xhelp2; + ihelp=48+(tof12_i)*3+tof31_i ; + c1 = k1[ihelp] ; + beta1[2] = c2*F/(ds-c1); +} +} + +//--------- S12 - S32 ------------------------------------------- + +if ((tof12_i>-1)&&(tof32_i>-1)) { + +dist = zin[1] - zin[5]; +F = 1./cos(theta); +c2 = (2.*0.01*dist)/(3.E08*50.E-12); + + ipmt[0] = 16+(tof12_i)*2; + ipmt[1] = 16+(tof12_i)*2+1; + ipmt[2] = 42+(tof32_i)*2; + ipmt[3] = 42+(tof32_i)*2+1; + + for (Int_t jj=0; jj< 4; jj++) time[jj] = tdca[(ipmt[jj])] ; + + if ((time[0]<4095)&&(time[1]<4095)&&(time[2]<4095)&&(time[3]<4095)) { + xhelp1 = time[0] + time[1] ; + xhelp2 = time[2] + time[3] ; + ds = xhelp1-xhelp2; + ihelp=66+(tof12_i)*3+tof32_i ; + c1 = k1[ihelp] ; + beta1[3] = c2*F/(ds-c1); +} +} + +//---------------------- calculate beta mean ----------------- + + sw=0.; + sxw=0.; + beta_mean_tof=100.; + + for (Int_t jj=0; jj<4;jj++){ + if ((beta1[jj]>0.1) && (beta1[jj]<1.5)) { + w_i=1./(0.13*0.13); + sxw=sxw + beta1[jj]*w_i ; + sw =sw + w_i ; + } + } + + if (sw>0) beta_mean_tof=sxw/sw; + +} // if tof11_i > -1 && ...... beta calculation + + + Float_t beta_help = beta_mean_tof ; // pow(beta_mean_tof,1.0) gave best results + +//----------------------- Select helium -------------------------- + + Int_t icount=0; + + for (jj=0; jj<24; jj++){ + a1 = adca[2*jj]*cos(theta); + a2 = adca[2*jj+1]*cos(theta); + + xhelp = 100000.; + if ((a1 < 3000) && (a2 < 3000)) xhelp = sqrt(a1*a2); // geometric mean +// if geometric mean multiplied by beta_help is inside helium limits, increase counter + if ((beta_mean_tof>0.6) && (beta_mean_tof<1.1) && + ((beta_help*xhelp)>A_l[jj]) && ((beta_help*xhelp) 3) iz=2; // if more than three paddles see helium, then set Z=2 + if (icount > 4) iz=2; + +//---------------------- Z=2 fill histograms ----------------------------- + + if (iz==2) { + + heevent++; + for (jj=0; jj<48; jj++) adche[jj]->Fill(adca[jj]); + + } // iz0==2 + + +//===================== end beta and helium part =========================== + + } // i < nevents + + float *X = new float[48]; float *means = new float[48]; float *entries = new float[48]; int *entriestdc = new int[48]; int *entriesadc = new int[48]; - + const char *saveas = format; - - // all range ADC/TDC - j = 0; - while (j < 12){ - k = 0; - z = 0; - - while(k < 4){ - if (k > 1) z = 2; - i = j*4 + k; - X[i] = i; - - tpTDCall = tdcHisto->ProjectionX(" ", i+1, i+1); - tpTDCall->SetAxisRange(0,4095); - entriestdc[i] = (Int_t)tpTDCall->Integral(); - - tpADCall = adcHisto->ProjectionX(" ", i+1, i+1); - tpADCall->SetAxisRange(0,4095); - entriesadc[i] = (Int_t)tpADCall->Integral(); - - k++; - } - - j++; - - } + int i=0; gStyle->SetStatW(0.4); gStyle->SetStatH(0.4); @@ -166,6 +464,7 @@ TCanvas *SCanvas = new TCanvas("SCanvas","SCanvas", 1280, 1024); SCanvas->Divide(4,2); + j = 0; while (j < 12){ k = 0; @@ -184,27 +483,31 @@ X[i] = i; SCanvas->cd(k+3+z); - tpTDC = tdcHisto->ProjectionX(" ", i+1, i+1); + htdc[i] = (TH1F*)hhtdc->At(i); + entriestdc[i] = (Int_t)htdc[i]->Integral(); sst.str(""); sst << "TDC - " << photoS[i].Data() << " (Nev < 4096 = " << entriestdc[i] << ")"; - tpTDC->SetTitle(sst.str().c_str()); - tpTDC->SetTitleSize(10); - tpTDC->SetAxisRange(690,1510); - tpTDC->DrawCopy(); - tpTDC->ComputeIntegral(); - entries[i] = tpTDC->Integral(); + htdc[i]->SetTitle(sst.str().c_str()); + htdc[i]->SetTitleSize(10); + htdc[i]->SetAxisRange(690,1510); + htdc[i]->DrawCopy(); + htdc[i]->ComputeIntegral(); + entries[i] = htdc[i]->Integral(); SCanvas->cd(k+1+z); - tpADC = adcHisto->ProjectionX(" ", i+1, i+1); + hadc[i] = (TH1F*)hhadc->At(i); + entriesadc[i] = (Int_t)hadc[i]->Integral(); sst.str(""); sst << "ADC - " << photoS[i].Data() << " (Nev < 4096 = " << entriesadc[i] << ")"; - tpADC->SetTitle(sst.str().c_str()); - tpADC->SetAxisRange(-10,710); - tpADC->DrawCopy(); - means[i] = tpADC->GetMean(); + hadc[i]->SetTitle(sst.str().c_str()); + hadc[i]->SetAxisRange(-10,710); + hadc[i]->DrawCopy(); + means[i] = hadc[i]->GetMean(); k++; } + + if ( !strcmp(saveas,"ps") ) { sst.str(""); sst << outDirectory.Data() << filename.Data() << "TOFScan.ps("; @@ -217,13 +520,15 @@ } j++; } + if (gROOT->IsBatch()) SCanvas->Close(); /* * This Canvas will represent a summary of the performances for TOF TDC/ADC channels */ - TCanvas *performanceCanvas = new TCanvas("performanceCanvas","performanceCanvas", 1280, 1024); - performanceCanvas->Divide(1,2); +// TCanvas *performanceCanvas = new TCanvas("performanceCanvas","performanceCanvas", 1280, 1024); + TCanvas *performanceCanvas = new TCanvas("performanceCanvas","performanceCanvas", 1024, 1280); + performanceCanvas->Divide(1,3); gStyle->SetTitleW(.9); @@ -232,7 +537,7 @@ sst.str(""); sst << "ADCMean" << " - Data in " << base.Data() << " - Nevents in the run = " << nevents; adcMeans->SetTitle(sst.str().c_str()); - adcMeans->SetFillColor(2); + adcMeans->SetFillColor(35); adcMeans->GetXaxis()->SetTitle("Photomultipliers"); adcMeans->GetXaxis()->CenterTitle(); adcMeans->GetXaxis()->SetLimits(-0.5, 47.5); @@ -245,14 +550,142 @@ sst.str(""); sst << "TDCEntries" << " - Data in " << base.Data() << " - Nevents in the run = " << nevents; tdcEntries->SetTitle(sst.str().c_str()); - tdcEntries->SetFillColor(2); + tdcEntries->SetFillColor(35); tdcEntries->GetXaxis()->SetTitle("Photomultipliers"); tdcEntries->GetXaxis()->CenterTitle(); tdcEntries->GetXaxis()->SetLimits(-0.5, 47.5); tdcEntries->GetYaxis()->SetTitle("TDCIntegral"); tdcEntries->GetYaxis()->CenterTitle(); tdcEntries->Draw("AB"); - + +//--------- new part PMT gain check ----------------------------- + + performanceCanvas->cd(3); + +Float_t xc[48],xmean1[48],xmeana[48]; +Float_t xmean_arr[12][48]; + +// xmean values from 2-3 april 2007 + +char date_info[]="Reference Data: apr-2007"; + +Float_t xmean[48] = { +491.609,509.241,400.786,530.122,699.674,555.747,521.04,486.363, +470.173,227.752,611.038,455.889,553.601,520.54,403.527,382.099, +349.697,365.113,447.653,377.667,517.815,572.932,338.501,436.681, +485.696,450.491,395.375,329.631,751.258,626.681,385.561,578.476, +374.454,356.733,641.888,562.767,582.849,521.748,527.043,505.89, +489.828,628.408,532.924,506.511,482.872,532.236,554.554,498.849 }; + +// new 01-oct-2007 +int channelmap[] = {0,7,3,6,2,8,1,5,3,7,3,6,1,7,2,10, + 10,10,10,5,0,7,0,5,0,6,1,5, + 2,8,3,8,2,6,1,8, + 11,9,11,11,9,11,4,4,4,9,9,4}; + + +int colormap[] = {46,2,29,4,5,6,7,8,9,11,28,34}; +//int colormap[] = {417,400,632,617,603,600,434,419,591,625,403,424}; + + +for (Int_t j=0; j<48; j++) xmeana[j]=0.; +for (Int_t j=0; j<24; j++) xmeana[2*j]=xmean[2*j]; + +for (Int_t i=0; i<12; i++) { +for (Int_t j=0; j<48; j++) { +xmean_arr[i][j]=0.; +} +} + +for (Int_t j=0; j<48; j++) { +Int_t ichan = channelmap[j]; +xmean_arr[ichan][j]=xmean[j]; +} + +// get results from ADC histogram +for (Int_t j=0; j<48; j++) { +xc[j]=j; +xmean1[j]=adche[j]->GetMean(); +} + + + gStyle->SetTitleW(.5); + gStyle->SetTitleH(.05); + +TH2F *hr = new TH2F("frame","2-Dim",2,-0.5,47.5,2,-300.,100.); +hr->SetStats(kFALSE); +hr->GetXaxis()->CenterTitle(); +hr->GetXaxis()->SetTitle("Photomultipliers"); +hr->GetYaxis()->CenterTitle(); +hr->GetYaxis()->SetTitle("Mean ADC Difference"); +hr->SetTitle("Difference between Reference and Actual Values"); +hr->Draw(); + +Int_t npoint=48; + +for (Int_t j=0; j<12; j++) { +for (Int_t i=0; i<48; i++) xmeana[i] = 0.; +for (Int_t i=0; i<48; i++) { +if (xmean_arr[j][i] != 0) xmeana[i] = xmean1[i] - xmean_arr[j][i]; + } + + +TGraph *graph1 = new TGraph(npoint,xc,xmeana); +graph1->SetFillColor(colormap[j]); +graph1->GetXaxis()->SetLimits(-0.5, 47.5); +graph1->Draw("BP"); +} + +Float_t tp[10]; +tp[0] = 15.5; +tp[1] = 27.5; +tp[2] = 31.5; +tp[3] = 35.5; +tp[4] = 41.5; + +for (Int_t ii=0; ii<5; ii++) { +TLine *l1=new TLine(tp[ii],-300,tp[ii],100); +l1->SetLineColor(38); +l1->Draw("same"); + } + +for (Int_t j=0; j<12; j++) { +sprintf(htitle, "HV_%d",j); +TText *text1 = new TText(0+j*4,80,htitle); +text1->SetTextColor(colormap[j]); +//text1->SetTextSize(0.03); +text1->SetTextSize(0.05); +text1->Draw(); +} + + +TText *text1 = new TText(0,-185,date_info); +text1->SetTextColor(kBlack); +text1->SetTextSize(0.023); +text1->Draw(); + + +sprintf(htitle, "Helium Events: %d",heevent); +TText *text2 = new TText(20,-185,htitle); +text2->SetTextColor(kBlack); +text2->SetTextSize(0.023); +text2->Draw(); + + +for (Int_t i=0; i<6; i++) { +for (Int_t j=0; j<8; j++) { +Int_t ihelp = i*8+j; +sprintf(htitle, "%d: %.0f/%.0f",(ihelp+1),xmean[ihelp],xmean1[ihelp]); +TText *text1 = new TText(0+j*6,-200-i*15,htitle); +text1->SetTextColor(kBlack); +text1->SetTextSize(0.023); +text1->Draw(); +} +} + +//-------- end new part ------------------------- + + //------print the ps if ( !strcmp(saveas,"ps") ) { @@ -270,6 +703,8 @@ performanceCanvas->Close(); } + + } int main(int argc, char* argv[]){ @@ -326,3 +761,4 @@ TofScan(argv[1], outDir, format); } +