Skip to content
Snippets Groups Projects
Commit 8011c768 authored by Simon Zhamkochyan's avatar Simon Zhamkochyan
Browse files

Added Thresholds, workong on adding ADC spectrums to the tree

parent f2a7e3c8
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
examples/compare.C 0 → 100644
+ 147
− 0
View file @ 8011c768
{
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();
*/
}
examples/hodtest.C
+ 2
− 2
View file @ 8011c768
......@@ -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" );
......
examples/output.def
+ 19
− 6
View file @ 8011c768
......@@ -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
src/THcScintillatorPlane.cxx
+ 3
− 3
View file @ 8011c768
......@@ -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);
......
src/THcShower.cxx
+ 1
− 1
View file @ 8011c768
......@@ -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;
......
src/THcShowerPlane.cxx
+ 108
− 7
View file @ 8011c768
......@@ -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);
}
......
src/THcShowerPlane.h
+ 15
− 0
View file @ 8011c768
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment