From 8011c768e72d0aa0d2dd4821418c3b035f48963c Mon Sep 17 00:00:00 2001 From: Simon Zhamkochyan <szh@jlab.org> Date: Tue, 7 Aug 2012 01:01:26 -0400 Subject: [PATCH] Added Thresholds, workong on adding ADC spectrums to the tree --- examples/compare.C | 147 +++++++++++++++++++++++++++++++++++ examples/hodtest.C | 4 +- examples/output.def | 25 ++++-- src/THcScintillatorPlane.cxx | 6 +- src/THcShower.cxx | 2 +- src/THcShowerPlane.cxx | 115 +++++++++++++++++++++++++-- src/THcShowerPlane.h | 15 ++++ 7 files changed, 295 insertions(+), 19 deletions(-) create mode 100644 examples/compare.C diff --git a/examples/compare.C b/examples/compare.C new file mode 100644 index 0000000..1664415 --- /dev/null +++ b/examples/compare.C @@ -0,0 +1,147 @@ +{ + TFile* f = new TFile("hodtest.root"); + + TCanvas *c1 = new TCanvas("c1", "Shower Hit Maps", 1000, 1200); + c1->Divide(2, 3); + + TH1F *h1[6]; + TH1F* h[6]; + + for(int l = 0; l < 6; l++){ + h1[l] = new TH1F("","",13,0.4,13.4); + } + Double_t mean = 0.; + Double_t hits; + Double_t raw[13] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13}; + Double_t ratio[77]; + Double_t err[77]; + + +//with peds and scalers +// double v[6][13] = { +// {679, 2016, 2439, 2446, 2369, 2179, 2045, 1953, 1723, 1618, 1465, 1241, 1061}, //hcalapos +// {0, 2125, 2521, 2420, 2364, 2006, 1935, 1911, 1769, 1536, 1450, 1242, 1090}, //hcalaneg +// {803, 2069, 2446, 2422, 2188, 1986, 2043, 1881, 1736, 1630, 1455, 1130, 1086}, //hcalbpos +// {788, 2034, 2372, 2217, 2141, 2002, 1956, 1938, 1700, 1605, 1442, 1219, 952}, //hcalbneg +// {915, 1964, 2276, 2066, 2122, 1880, 1969, 1854, 1576, 1500, 1403, 1189, 875}, //hcalcpos +// {7187, 1748, 1890, 1852, 1786, 1573, 1525, 1510, 1376, 1432, 1191, 1033, 809} //hcaldpos +// }; + + double v[6][13] = { + {747, 2212, 2688, 2698, 2602, 2400, 2228, 2131, 1910, 1763, 1605, 1382, 1179}, //hcalapos + {0, 2334, 2784, 2673, 2599, 2215, 2107, 2079, 1958, 1672, 1584, 1384, 1211}, //hcalaneg + {882, 2285, 2709, 2650, 2405, 2184, 2230, 2062, 1921, 1789, 1589, 1252, 1204}, //hcalbpos + {866, 2244, 2632, 2430, 2354, 2207, 2127, 2127, 1888, 1757, 1578, 1345, 1054}, //hcalbneg + {1006, 2171, 2529, 2262, 2302, 2070, 2149, 2024, 1732, 1634, 1543, 1323, 957}, //hcalcpos + {7913, 1918, 2084, 2027, 1947, 1739, 1668, 1646, 1521, 1572, 1302, 1136, 887} //hcaldpos + }; + + + // //Without thresholds + // double v[6][13] = { + // {679, 2296, 2439, 2446, 2369, 2179, 2045, 1953, 1723, 1618, 1465, 1241, 1061}; + // {0, 2125, 2521, 2420, 2364, 2006, 1935, 1911, 1769, 1536, 1450, 1242, 1090}; + // {803, 2069, 2446, 2422, 2188, 1986, 2043, 1881, 1736, 1630, 1455, 1130, 1086}; + // {788, 2034, 2372, 2217, 2141, 2002, 1956, 1938, 1700, 1605, 1442, 1219, 952}; + // {915, 1964, 2276, 2066, 2122, 1880, 1969, 1854, 1576, 1500, 1403, 1189, 875}; + // {12633, 1748, 1890, 1852, 1786, 1573, 1525, 1510, 1376, 1432, 1191, 1033, 809}; + + h[0] = chposadc1; + h[1] = chnegadc1; + h[2] = chposadc2; + h[3] = chnegadc2; + h[4] = chposadc3; + h[5] = chposadc4; + + + for (int j = 0; j < 6; j++){ + c1->cd(j+1); + h[j]->SetFillColor(kGreen); + h[j]->SetFillStyle(3345); + h[j]->Draw(); + + for(int k = 0; k < 13; k++){ + for(int i = 0; i < v[j][k]; i++){ + hits = raw[k]-0.1; + h1[j]->Fill(hits); + } + } + h1[j]->SetFillColor(kBlue); + h1[j]->SetFillStyle(3354); + h1[j]->Draw("same"); + } + +for (int i = 1; i < 14; i++){ +cout << "hitsC = " << h[0].GetBinContent(i) << endl; +} + + +/* +int n = 0; +Double_t errC; +Double_t errE; +Double_t hitsC; +Double_t hitsE; +Double_t hitsCmean = 1876.81; +Double_t hitsEmean = 1783.7; +Double_t sum1 = 0.; //chislitel +Double_t sum2 = 0.; +Double_t sum3 = 0.; +Double_t chan[77]; +Double_t r; + + for (Int_t l = 0; l < 6; l++){ +// cout << "ADC" << l+1 << ": " << endl; + for (Int_t i = 1; i < h[l].GetNbinsX() +1; i++){ + hitsC = h[l].GetBinContent(i); + hitsE = h1[l].GetBinContent(i); + if (hitsC > 0 && hitsE > 0){ + errC = 1./sqrt(hitsC); + errE = 1./sqrt(hitsE); + ratio[n] = hitsC/hitsE; +//cout << "errC = " << errC << ", errE = " << errE << "corr = " << errC*errE << endl; + err[n] = sqrt((errC**2 + errE**2 - 2*errC*errE))*ratio[n]; + mean = mean + ratio[n]; + + + + +sum1 = sum1 + (hitsC - hitsCmean)*(hitsE - hitsEmean); +sum2 = sum2 + (hitsC - hitsCmean)**2; +sum3 = sum3 + (hitsE - hitsEmean)**2; +// hitsCmean = hitsCmean + hitsC; +// hitsEmean = hitsEmean + hitsE; + chan[n] = n*1.; + n++; + } + //cout << "Hits ratio C++/ENGINE = " << err[l][i] << endl; + } + } +mean = mean/n; +//hitsCmean = hitsCmean/n; +//hitsEmean = hitsEmean/n; +r = sum1/sqrt(sum2*sum3); +cout << "r = " << r << endl; + +c1->Update(); +TCanvas *c2 = new TCanvas("c2", "RATIO", 1200, 800); +c2->cd(); +// create the TGraphErrors and draw it +gr = new TGraphErrors(n,chan,ratio,0,err); +//gr->SetLineColor(8); +gr->SetMarkerStyle(20); + TAxis *xaxis = gr->GetXaxis(); + TAxis *yaxis = gr->GetYaxis(); + //xaxis->SetLimits(10, 60); + xaxis->SetTitle("Channel"); + yaxis->SetTitle("Hits ratio (C++/Engine)"); +yaxis->SetTitleOffset(1.2); + +gr->Draw("AP"); + + TF1 *f1 = new TF1("f1","1.0527",0.,77.); + f1->Draw("same"); + +c2->Update(); +*/ +} diff --git a/examples/hodtest.C b/examples/hodtest.C index 2ee29bb..240160b 100644 --- a/examples/hodtest.C +++ b/examples/hodtest.C @@ -39,8 +39,8 @@ gHcDetectorMap->Load("jan03.map"); // We just set up one, but this could be many. THaRun* run = new THaRun( "daq04_50017.log.0" ); //THaRun* run = new THaRun( "daq03_47851.log.0" ); - run->SetEventRange(2000,100000); - + run->SetEventRange(1054,100000); + // Define the analysis parameters analyzer->SetEvent( event ); analyzer->SetOutFile( "hodtest.root" ); diff --git a/examples/output.def b/examples/output.def index cf1c55c..194c7ce 100644 --- a/examples/output.def +++ b/examples/output.def @@ -26,10 +26,23 @@ TH1F hposadc4 'HMS s2y+ ADC hits' H.hod.2y.posadchits 10 0.5 10.5 TH1F hnegadc4 'HMS s2y- ADC hits' H.hod.2y.negadchits 10 0.5 10.5 # ADC hits per Calorimeter layer -TH1F chposadc1 'HMS Cal 1z+ ADC hits' H.Cal.1z.posadchits 13 0.5 13.5 -TH1F chnegadc1 'HMS Cal 1z- ADC hits' H.Cal.1z.negadchits 13 0.5 13.5 -TH1F chposadc2 'HMS Cal 2z+ ADC hits' H.Cal.2z.posadchits 13 0.5 13.5 -TH1F chnegadc2 'HMS Cal 2z- ADC hits' H.Cal.2z.negadchits 13 0.5 13.5 -TH1F chposadc3 'HMS Cal 3z+ ADC hits' H.Cal.3z.posadchits 13 0.5 13.5 -TH1F chposadc4 'HMS Cal 4z+ ADC hits' H.Cal.4z.posadchits 13 0.5 13.5 +TH1F chposadc1 'HMS Cal A+ ADC hits' H.Cal.1z.posadchits 13 0.5 13.5 +TH1F chnegadc1 'HMS Cal A- ADC hits' H.Cal.1z.negadchits 13 0.5 13.5 +TH1F chposadc2 'HMS Cal B+ ADC hits' H.Cal.2z.posadchits 13 0.5 13.5 +TH1F chnegadc2 'HMS Cal B- ADC hits' H.Cal.2z.negadchits 13 0.5 13.5 +TH1F chposadc3 'HMS Cal C+ ADC hits' H.Cal.3z.posadchits 13 0.5 13.5 +TH1F chposadc4 'HMS Cal D+ ADC hits' H.Cal.4z.posadchits 13 0.5 13.5 +TH1F calposadc1 'HMS Cal ADC1' H.Cal.1z.posadc1 150 -50 400 +# TH1F calposadc2 'HMS Cal ADC2' H.Cal.2z.posadc1 150 -50 400 +# TH1F chposadc 'HMS Cal ADC3' H.Cal.1z.posadc3 150 -50 400 +# TH1F chposadc 'HMS Cal ADC4' H.Cal.1z.posadc4 150 -50 400 +# TH1F chposadc 'HMS Cal ADC5' H.Cal.1z.posadc5 150 -50 400 +# TH1F chposadc 'HMS Cal ADC6' H.Cal.1z.posadc6 150 -50 400 +# TH1F chposadc 'HMS Cal ADC7' H.Cal.1z.posadc7 150 -50 400 +# TH1F chposadc 'HMS Cal ADC8' H.Cal.1z.posadc8 150 -50 400 +# TH1F chposadc 'HMS Cal ADC9' H.Cal.1z.posadc9 150 -50 400 +# TH1F chposadc 'HMS Cal ADC10' H.Cal.1z.posadc10 150 -50 400 +# TH1F chposadc 'HMS Cal ADC11' H.Cal.1z.posadc11 150 -50 400 +# TH1F chposadc 'HMS Cal ADC12' H.Cal.1z.posadc12 150 -50 400 +# TH1F chposadc 'HMS Cal ADC13' H.Cal.1z.posadc13 150 -50 400 \ No newline at end of file diff --git a/src/THcScintillatorPlane.cxx b/src/THcScintillatorPlane.cxx index 4fe536e..7f966c2 100644 --- a/src/THcScintillatorPlane.cxx +++ b/src/THcScintillatorPlane.cxx @@ -231,7 +231,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) THcSignalHit *sighit = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++); #else - TObject* obj = (*fPosTDCHits)[nNegTDCHits++]; + TObject* obj = (*fNegTDCHits)[nNegTDCHits++]; R__ASSERT( obj ); if(!obj->TestBit (TObject::kNotDeleted)) fNegTDCHitsClass->New(obj); @@ -244,7 +244,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) #if ROOT_VERSION_CODE >= ROOT_VERSION(5,32,0) THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++); #else - TObject* obj = (*fPosTDCHits)[nPosADCHits++]; + TObject* obj = (*fPosADCHits)[nPosADCHits++]; R__ASSERT( obj ); if(!obj->TestBit (TObject::kNotDeleted)) fPosADCHitsClass->New(obj); @@ -257,7 +257,7 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) #if ROOT_VERSION_CODE >= ROOT_VERSION(5,32,0) THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++); #else - TObject* obj = (*fPosTDCHits)[nNegADCHits++]; + TObject* obj = (*fNegADCHits)[nNegADCHits++]; R__ASSERT( obj ); if(!obj->TestBit (TObject::kNotDeleted)) fNegADCHitsClass->New(obj); diff --git a/src/THcShower.cxx b/src/THcShower.cxx index fe58c89..a551db5 100644 --- a/src/THcShower.cxx +++ b/src/THcShower.cxx @@ -90,7 +90,7 @@ THaAnalysisObject::EStatus THcShower::Init( const TDatime& date ) { static const char* const here = "Init()"; - cout << "THcHodoscope::Init " << GetName() << endl; + cout << "THcShower::Init " << GetName() << endl; if( THaNonTrackingDetector::Init( date ) ) return fStatus; diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx index 0fc4e4e..42dbeda 100644 --- a/src/THcShowerPlane.cxx +++ b/src/THcShowerPlane.cxx @@ -20,6 +20,7 @@ #include <cstdlib> #include <iostream> +#include <fstream> using namespace std; ClassImp(THcShowerPlane) @@ -32,10 +33,14 @@ THcShowerPlane::THcShowerPlane( const char* name, : THaSubDetector(name,description,parent) { // Normal constructor with name and description - fPosADCHits = new TClonesArray("THcSignalHit",16); - fNegADCHits = new TClonesArray("THcSignalHit",16); + fPosADCHits = new TClonesArray("THcSignalHit",13); + fNegADCHits = new TClonesArray("THcSignalHit",13); fPosADCHitsClass = fPosADCHits->GetClass(); fNegADCHitsClass = fNegADCHits->GetClass(); + + fPosADC1 = new TClonesArray("THcSignalHit",13); + fPosADCHitsClass = fPosADC1->GetClass(); + fLayerNum = layernum; } @@ -45,6 +50,7 @@ THcShowerPlane::~THcShowerPlane() // Destructor delete fPosADCHits; delete fNegADCHits; + delete fPosADC1; } THaAnalysisObject::EStatus THcShowerPlane::Init( const TDatime& date ) @@ -116,15 +122,79 @@ Int_t THcShowerPlane::DefineVariables( EMode mode ) if( mode == kDefine && fIsSetup ) return kOK; fIsSetup = ( mode == kDefine ); +/* +char spos1[256]; +char spos2[256]; +char spos3[256]; + + // Register variables in global list + RVarDef vars[15]; + vars[0] = "posadchits", "List of Positive ADC hits", + "fPosADCHits.THcSignalHit.GetPaddleNumber()"; + + vars[1] = "negadchits", "List of Negative ADC hits", + "fNegADCHits.THcSignalHit.GetPaddleNumber()"; + + for(Int_t i=2;i<15;i++){ + sprintf(spos1,"posadc1%d",i-1); + sprintf(spos2,"ADC%d hits",i-1); + sprintf(spos3,"fPosADC1[%d].THcSignalHit.Get",i-1); + vars[i] = spos1,spos2,spos3; + } + +*/ +fA = new Float_t[13]; + +char spos1[256]; +char spos2[256]; +char spos3[256]; + +sprintf(spos1,"posadc%d",1); +sprintf(spos2,"ADC%d hits",1); +sprintf(spos3,"fPosADC%d.THcSignalHit.GetPaddleNumber()",1); +Double_t *fa = new Double_t[4]; +fa[0] = 0.0; +fa[1] = 1.0; +fa[2] = 2.0; +fa[3] = 3.0; + +RVarDef vars[4]; + +vars[0].name = "posadchits"; +vars[0].desc = "List of Positive ADC hits"; +vars[0].def = "fPosADCHits.THcSignalHit.GetPaddleNumber()"; + +vars[1].name = "negadchits"; +vars[1].desc = "List of Negative ADC hits"; +vars[1].def = "fNegADCHits.THcSignalHit.GetPaddleNumber()"; + +vars[2].name = "posadc1"; +vars[2].desc = "ADC1 hits"; +//vars[2].def = "fA"; +vars[2].def = "fPosADC1.THcSignalHit.GetPaddleNumber()"; + +// vars[2].name = spos1; +// vars[2].desc = spos2; +// vars[2].def = spos3; + + +vars[3].name = NULL; +vars[3].desc = NULL; +vars[3].def = NULL; + + +/* // Register variables in global list RVarDef vars[] = { {"posadchits", "List of Positive ADC hits", "fPosADCHits.THcSignalHit.GetPaddleNumber()"}, {"negadchits", "List of Negative ADC hits", "fNegADCHits.THcSignalHit.GetPaddleNumber()"}, + {"posadc1", "ADC1 hits", + "fPosADC1.THcSignalHit.GetPaddleNumber()"}, { 0 } }; - +*/ return DefineVarsFromList( vars, mode ); } @@ -135,6 +205,7 @@ void THcShowerPlane::Clear( Option_t* ) // Clears the hit lists fPosADCHits->Clear(); fNegADCHits->Clear(); + fPosADC1->Clear(); } //_____________________________________________________________________________ @@ -169,7 +240,27 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) Int_t nNegADCHits=0; fPosADCHits->Clear(); fNegADCHits->Clear(); + fPosADC1->Clear(); + +//-------Thresholds values are taken from ENGINE +double hcal_new_threshold_pos[52] = { +// 485.7,505.7,407.8,513.6,547.2,343.0,437.5,420.2,444.6,354.7,498.4,251.9,565.3, +// 430.5,485.1,286.7,388.4,425.9,452.1,486.8,456.3,285.9,381.8,250.9,261.6,424.7, +// 361.4,378.9,431.3,224.0,488.2,352.9,436.3,379.9,343.9,432.8,462.8,252.3,427.0, +// 399.9,296.9,378.6,367.4,514.9,564.8,561.5,475.8,390.1,433.9,334.4,407.4,516.2}; +470.7+12.6, 490.7+24.1, 392.8+12.1, 498.6+10.9, 532.2+14.1, 328.0+14.5, 422.5+14.0, 405.2+12.5, 429.6+11.4, 339.7+12.9, 483.4+14.8, 236.9+11.5, 550.3+10.0, +415.5+12.5, 470.1+16.3, 271.7+11.7, 373.4+11.5, 410.9+10.7, 437.1+10.3, 471.8+10.0, 441.3+10.0, 270.9+11.0, 366.8+10.0, 235.9+12.5, 246.6+10.0, 409.7+11.6, +346.4+13.2, 363.9+11.1, 416.3+10.1, 209.0+10.0, 473.2+12.5, 337.9+14.4, 421.3+10.0, 364.9+10.1, 328.9+17.9, 417.8+11.6, 447.8+16.5, 237.3+11.5, 412.0+11.4, +349.9+50.0, 281.9+11.2, 363.6+18.1, 352.4+14.7, 499.9+14.7, 549.8+15.3, 546.5+14.6, 460.8+13.6, 375.1+13.4, 418.9+11.0, 319.4+14.1, 392.4+12.2, 501.2+16.0}; + + +double hcal_new_threshold_neg[52] = { +520.8+10.0, 472.0+15.9, 450.6+15.1, 451.5+12.9, 523.0+14.0, 553.3+15.9, 573.0+17.7, 494.9+14.3, 493.8+16.2, 487.8+12.2, 436.5+16.1, 438.0+10.0, 573.1+13.5, +527.4+11.8, 465.0+13.7, 460.1+13.6, 390.8+12.7, 552.3+16.7, 623.3+13.1, 549.3+13.3, 632.9+10.0, 479.7+12.1, 481.9+10.3, 412.3+15.0, 410.4+11.8, 635.1+12.6}; + +CalADC1File = fopen("adc1_new.dat", "a"); +//fprintf(CalADC1File, "%d\n", 1); Int_t nrawhits = rawhits->GetLast()+1; Int_t ihit = nexthit; @@ -179,11 +270,18 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) break; } +if(hit->fCounter == 1){ +THcSignalHit *sighit1 = (THcSignalHit*) fPosADC1->ConstructedAt(nPosADCHits++); +//THcSignalHit *sighit1 = (THcSignalHit*) fA[1]->ConstructedAt(nPosADCHits++); +sighit1->Set(hit->fADC_pos - 470.7,1); +//fprintf(CalADC1File, "%d\n", hit->fADC_pos); +} -if(hit->fADC_pos > 0) { +double thresh_pos = hcal_new_threshold_pos[hit->fCounter + 13*(hit->fPlane -1) -1]; +if(hit->fADC_pos > thresh_pos) { #if ROOT_VERSION_CODE >= ROOT_VERSION(5,32,0) THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++); - sighit->Set(hit->fCounter, hit->fADC_pos); + #else TObject* obj = (*fPosADCHits)[nPosADCHits++]; R__ASSERT( obj ); @@ -191,15 +289,17 @@ if(!obj->TestBit (TObject::kNotDeleted)) fPosADCHitsClass->New(obj); THcSignalHit *sighit = (THcSignalHit*)obj; #endif + sighit->Set(hit->fCounter, hit->fADC_pos); } -if(hit->fADC_neg > 0) { +double thresh_neg = hcal_new_threshold_neg[hit->fCounter + 13*(hit->fPlane -1) -1]; +if(hit->fADC_neg > thresh_neg) { #if ROOT_VERSION_CODE >= ROOT_VERSION(5,32,0) THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++); sighit->Set(hit->fCounter, hit->fADC_neg); #else - TObject* obj = (*fPosADCHits)[nNegADCHits++]; + TObject* obj = (*fNegADCHits)[nNegADCHits++]; R__ASSERT( obj ); if(!obj->TestBit (TObject::kNotDeleted)) fNegADCHitsClass->New(obj); @@ -210,6 +310,7 @@ if(!obj->TestBit (TObject::kNotDeleted)) ihit++; } +fclose(CalADC1File); return(ihit); } diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h index ec0d7eb..4e44d1c 100644 --- a/src/THcShowerPlane.h +++ b/src/THcShowerPlane.h @@ -15,6 +15,10 @@ #include "THaSubDetector.h" #include "TClonesArray.h" +#include <iostream> + +#include <fstream> + class THaEvData; class THaSignalHit; @@ -42,11 +46,22 @@ class THcShowerPlane : public THaSubDetector { protected: + +// TClonesArray* fPosADC1[13]; + + + Float_t* fA; // [fNelem] Array of ADC amplitudes of blocks +TClonesArray* fPosADC1; +TClonesArray* fPosADC[13]; + TClonesArray* fPosADCHits; TClonesArray* fNegADCHits; TClass* fPosADCHitsClass; TClass* fNegADCHitsClass; + TClass* fPosADC1Class; + + FILE* CalADC1File; Int_t fLayerNum; -- GitLab