diff --git a/examples/PARAM/52949/hcana.param b/examples/PARAM/52949/hcana.param
index 61d60c67416cd2b92a6e69416e20b1281f30689b..d49faaa41e6d7cb2de2eae66d474da35c8f4c7d1 100644
--- a/examples/PARAM/52949/hcana.param
+++ b/examples/PARAM/52949/hcana.param
@@ -7,6 +7,8 @@ hhodo_num_planes = 4
hcal_num_layers = 4
+hcal_fv_delta = 5.
+
haero_num_pairs = 8
# Names of planes so that parameter names can be constructed
@@ -15,3 +17,6 @@ hdc_plane_names = "1x1 1y1 1u1 1v1 1y2 1x2 2x1 2y1 2u1 2v1 2y2 2x2"
hcal_layer_names = "1pr 2ta 3ta 4ta"
hhodo_plane_names = "1x 1y 2x 2y"
+
+hpcentral=0.811
+
diff --git a/examples/comptrackedep.C b/examples/comptrackedep.C
new file mode 100644
index 0000000000000000000000000000000000000000..d9ae5997b7b19eaa913a6db16c11997dba74331e
--- /dev/null
+++ b/examples/comptrackedep.C
@@ -0,0 +1,57 @@
+void comptrackedep(Int_t run)
+{
+ TFile* f = new TFile(Form("hodtest_%d.root",run));
+ cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+ TH1F* h = tre;
+
+ TFile* f1 = new TFile(Form("%d_hbk.root",run));
+ cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+ TH1F* h1;
+ switch (run) {
+ case 50017 :
+ // h1 = h212; //A+
+ break;
+ default :
+ h1 = h420; //edep
+ }
+
+ TCanvas *c1 = new TCanvas("c1", "Shower 1 track cluster Edep", 1000, 667);
+
+ // gPad->SetLogy();
+
+ h1->SetFillColor(kGreen);
+ h1->SetLineColor(kGreen);
+ h1->Draw();
+
+ h->SetLineColor(kBlue);
+ h->SetFillStyle(0);
+ h->SetLineWidth(2);
+ h->Draw("same");
+
+ TLatex l;
+ l.SetTextSize(0.04);
+ Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+ Float_t xmin = h1->GetXaxis()->GetXmin();
+ Float_t xmax = h1->GetXaxis()->GetXmax();
+ Float_t xt = xmin + 0.67*(xmax-xmin);
+
+ l.SetTextColor(kGreen);
+ l.DrawLatex(xt,0.65*maxy,"Engine");
+ l.SetTextColor(kBlue);
+ l.DrawLatex(xt,0.75*maxy,"hcana");
+
+ // Difference between the histograms.
+
+ TCanvas *c2 = new TCanvas("c2", "Edep differences", 1000, 667);
+
+ TH1F* dif = h->Clone();
+
+ dif->Add(h,h1,1.,-1.);
+
+ dif->SetTitle("Edep Difference");
+ dif->SetFillColor(kRed);
+ dif->SetLineColor(kRed);
+ dif->SetLineWidth(1);
+ dif->SetFillStyle(1111);
+ dif->Draw();
+}
diff --git a/examples/comptrackepr.C b/examples/comptrackepr.C
new file mode 100644
index 0000000000000000000000000000000000000000..a02a6e1d5cf88158e35d62694ab39e6bb76ae179
--- /dev/null
+++ b/examples/comptrackepr.C
@@ -0,0 +1,57 @@
+void comptrackepr(Int_t run)
+{
+ TFile* f = new TFile(Form("hodtest_%d.root",run));
+ cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+ TH1F* h = trepr;
+
+ TFile* f1 = new TFile(Form("%d_hbk.root",run));
+ cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+ TH1F* h1;
+ switch (run) {
+ case 50017 :
+ // h1 = h212; //A+
+ break;
+ default :
+ h1 = h421; //epr
+ }
+
+ TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster Edep", 1000, 667);
+
+ // gPad->SetLogy();
+
+ h1->SetFillColor(kGreen);
+ h1->SetLineColor(kGreen);
+ h1->Draw();
+
+ h->SetLineColor(kBlue);
+ h->SetFillStyle(0);
+ h->SetLineWidth(2);
+ h->Draw("same");
+
+ TLatex l;
+ l.SetTextSize(0.04);
+ Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+ Float_t xmin = h1->GetXaxis()->GetXmin();
+ Float_t xmax = h1->GetXaxis()->GetXmax();
+ Float_t xt = xmin + 0.75*(xmax-xmin);
+
+ l.SetTextColor(kGreen);
+ l.DrawLatex(xt,0.65*maxy,"Engine");
+ l.SetTextColor(kBlue);
+ l.DrawLatex(xt,0.75*maxy,"hcana");
+
+ // Difference between the histograms.
+
+ TCanvas *c2 = new TCanvas("c2", "Epr differences", 1000, 667);
+
+ TH1F* dif = h->Clone();
+
+ dif->Add(h,h1,1.,-1.);
+
+ dif->SetTitle("Epr Difference");
+ dif->SetFillColor(kRed);
+ dif->SetLineColor(kRed);
+ dif->SetLineWidth(1);
+ dif->SetFillStyle(1111);
+ dif->Draw();
+}
diff --git a/examples/comptrackx.C b/examples/comptrackx.C
new file mode 100644
index 0000000000000000000000000000000000000000..c67aac2cae708bde04427bdc0a1537c2ab8858a3
--- /dev/null
+++ b/examples/comptrackx.C
@@ -0,0 +1,57 @@
+void comptrackx(Int_t run)
+{
+ TFile* f = new TFile(Form("hodtest_%d.root",run));
+ cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+ TH1F* h = trx;
+
+ TFile* f1 = new TFile(Form("%d_hbk.root",run));
+ cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+ TH1F* h1;
+ switch (run) {
+ case 50017 :
+ // h1 = h212; //A+
+ break;
+ default :
+ h1 = h417; //X
+ }
+
+ TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster X", 1000, 667);
+
+ // gPad->SetLogy();
+
+ h1->SetFillColor(kGreen);
+ h1->SetLineColor(kGreen);
+ h1->Draw();
+
+ h->SetLineColor(kBlue);
+ h->SetFillStyle(0);
+ h->SetLineWidth(2);
+ h->Draw("same");
+
+ TLatex l;
+ l.SetTextSize(0.04);
+ Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+ Float_t xmin = h1->GetXaxis()->GetXmin();
+ Float_t xmax = h1->GetXaxis()->GetXmax();
+ Float_t xt = xmin + 0.67*(xmax-xmin);
+
+ l.SetTextColor(kGreen);
+ l.DrawLatex(xt,0.65*maxy,"Engine");
+ l.SetTextColor(kBlue);
+ l.DrawLatex(xt,0.75*maxy,"hcana");
+
+ // Difference between the histograms.
+
+ TCanvas *c2 = new TCanvas("c2", "Track X differences", 1000, 667);
+
+ TH1F* dif = h->Clone();
+
+ dif->Add(h,h1,1.,-1.);
+
+ dif->SetTitle("Hcal track X Difference");
+ dif->SetFillColor(kRed);
+ dif->SetLineColor(kRed);
+ dif->SetLineWidth(1);
+ dif->SetFillStyle(1111);
+ dif->Draw();
+}
diff --git a/examples/comptrackxy.C b/examples/comptrackxy.C
new file mode 100644
index 0000000000000000000000000000000000000000..b0c4b910235d62c3c6494858bad0b5a4d4184607
--- /dev/null
+++ b/examples/comptrackxy.C
@@ -0,0 +1,59 @@
+void comptrackxy(Int_t run)
+{
+ TFile* f = new TFile(Form("hodtest_%d.root",run));
+ cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+ TH2F* h = trx_vs_try;
+
+ TFile* f1 = new TFile(Form("%d_hbk.root",run));
+ cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+ TH2F* h1;
+ switch (run) {
+ case 50017 :
+ // h1 = h212; //A+
+ break;
+ default :
+ h1 = h419; //Y vs Y
+ }
+
+ TCanvas *c1 = new TCanvas("c1", "Shower Largest cluster X", 1000, 667);
+
+ // gPad->SetLogy();
+
+ h1->SetMarkerColor(kGreen);
+ // h1->SetLineColor(kGreen);
+ h1->Draw();
+
+ h->SetMarkerColor(kBlue);
+ // h->SetFillStyle(0);
+ // h->SetLineWidth(2);
+ h->Draw("same");
+
+ /*
+ TLatex l;
+ l.SetTextSize(0.04);
+ Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+ Float_t xmin = h1->GetXaxis()->GetXmin();
+ Float_t xmax = h1->GetXaxis()->GetXmax();
+ Float_t xt = xmin + 0.67*(xmax-xmin);
+
+ l.SetTextColor(kGreen);
+ l.DrawLatex(xt,0.095*maxy,"Engine");
+ l.SetTextColor(kBlue);
+ l.DrawLatex(xt,0.045*maxy,"hcana");
+ */
+
+ // Difference between the histograms.
+
+ TCanvas *c2 = new TCanvas("c2", "X differences", 1000, 667);
+
+ TH1F* dif = h->Clone();
+
+ dif->Add(h,h1,1.,-1.);
+
+ // dif->SetTitle("X Difference");
+ // dif->SetFillColor(kRed);
+ // dif->SetLineColor(kRed);
+ // dif->SetLineWidth(1);
+ // dif->SetFillStyle(1111);
+ dif->Draw("LEGO2");
+}
diff --git a/examples/comptracky.C b/examples/comptracky.C
new file mode 100644
index 0000000000000000000000000000000000000000..554a518b2f99d7e1d915bd438e9e57eb98131ba3
--- /dev/null
+++ b/examples/comptracky.C
@@ -0,0 +1,57 @@
+void comptracky(Int_t run)
+{
+ TFile* f = new TFile(Form("hodtest_%d.root",run));
+ cout << "hcana root file " << Form("hodtest_%d.root",run) << endl;
+ TH1F* h = try;
+
+ TFile* f1 = new TFile(Form("%d_hbk.root",run));
+ cout << "Engine root file " << Form("%d_hbk.root",run) << endl;
+ TH1F* h1;
+ switch (run) {
+ case 50017 :
+ // h1 = h212; //A+
+ break;
+ default :
+ h1 = h418; //Y
+ }
+
+ TCanvas *c1 = new TCanvas("c1", "Shower track Y", 1000, 667);
+
+ // gPad->SetLogy();
+
+ h1->SetFillColor(kGreen);
+ h1->SetLineColor(kGreen);
+ h1->Draw();
+
+ h->SetLineColor(kBlue);
+ h->SetFillStyle(0);
+ h->SetLineWidth(2);
+ h->Draw("same");
+
+ TLatex l;
+ l.SetTextSize(0.04);
+ Float_t maxy = h1->GetBinContent(h1->GetMaximumBin());
+ Float_t xmin = h1->GetXaxis()->GetXmin();
+ Float_t xmax = h1->GetXaxis()->GetXmax();
+ Float_t xt = xmin + 0.67*(xmax-xmin);
+
+ l.SetTextColor(kGreen);
+ l.DrawLatex(xt,0.65*maxy,"Engine");
+ l.SetTextColor(kBlue);
+ l.DrawLatex(xt,0.755*maxy,"hcana");
+
+ // Difference between the histograms.
+
+ TCanvas *c2 = new TCanvas("c2", "Hcal Y differences", 1000, 667);
+
+ TH1F* dif = h->Clone();
+
+ dif->Add(h,h1,1.,-1.);
+
+ dif->SetTitle("Y Difference");
+ dif->SetFillColor(kRed);
+ dif->SetLineColor(kRed);
+ dif->SetLineWidth(1);
+ dif->SetFillStyle(1111);
+ dif->Draw();
+}
diff --git a/examples/output_52949.def b/examples/output_52949.def
index a59a141addf4ece13baf3eb1ca330be31cbdd599..283b044a4bcd26a8b73cd4b04ed259ef8c1ffbe9 100644
--- a/examples/output_52949.def
+++ b/examples/output_52949.def
@@ -63,6 +63,17 @@ TH1F eprmax 'HMS PrSh Max Clust. Edep' H.cal.eprmax 100 -0.1 0.811
# X coordinate the largest cluster.
TH1F xmax 'HMS Cal Max Clust. X' H.cal.xmax 150 -75. 75.
+# Track coordinates at the front of calorimeter
+TH2F trx_vs_try 'HMS Track X vs Y at calor.' H.cal.try H.cal.trx 80 -40. 40. 150 -75. 75. H.dc.ntrack==1
+TH1F trx 'HMS Track X at calorimeter' H.cal.trx 150 -75. 75. H.dc.ntrack==1
+TH1F try 'HMS Track Y at calorimeter' H.cal.try 80 -40. 40. H.dc.ntrack==1
+
+# Energy of the cluster associated to the main (currently first) track
+TH1F tre 'HMS Cal Track Clust. Edep' H.cal.tre 200 -0.1 1.622 H.dc.ntrack==1
+
+# Preshower Energy of the cluster associated to the main (currently 1st) track
+TH1F trepr 'HMS PrSh Track Clust. Edep' H.cal.trepr 100 -0.1 0.811 H.dc.ntrack==1
+
#Calorimeter ADC channels
TH1F hcaladc_A1p 'HMS Cal ADC A1p - PED' H.cal.1pr.apos_p[0] 150 50 500
TH1F hcaladc_A2p 'HMS Cal ADC A2p - PED' H.cal.1pr.apos_p[1] 150 50 500
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index 4c8aff026a5f3a572609da67301f76084ed58029..f57f2e40a6739d8e7c0a35c112e80c406082171d 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -74,7 +74,9 @@ void THcShower::Setup(const char* name, const char* description)
vector<string> layer_names = vsplit(layernamelist);
if(layer_names.size() != (UInt_t) fNLayers) {
- cout << "ERROR: Number of layers " << fNLayers << " doesn't agree with number of layer names " << layer_names.size() << endl;
+ cout << "ERROR: Number of layers " << fNLayers
+ << " doesn't agree with number of layer names "
+ << layer_names.size() << endl;
// Should quit. Is there an official way to quit?
}
@@ -165,17 +167,25 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
DBRequest list[]={
{"cal_num_neg_columns", &fNegCols, kInt},
{"cal_slop", &fSlop, kDouble},
- {"cal_fv_test", &fvTest, kDouble},
+ {"cal_fv_test", &fvTest, kInt},
+ {"cal_fv_delta", &fvDelta, kDouble},
+ {"dbg_decoded_cal", &fdbg_decoded_cal, kInt},
+ {"dbg_sparsified_cal", &fdbg_sparsified_cal, kInt},
{"dbg_clusters_cal", &fdbg_clusters_cal, kInt},
+ {"dbg_tracks_cal", &fdbg_tracks_cal, kInt},
{0}
};
gHcParms->LoadParmValues((DBRequest*)&list, prefix);
}
- cout << "Number of neg. columns = " << fNegCols << endl;
- cout << "Slop parameter = " << fSlop << endl;
- cout << "Fiducial volum test flag = " << fvTest << endl;
- cout << "Cluster debug flag = " << fdbg_clusters_cal << endl;
+ cout << "Number of neg. columns = " << fNegCols << endl;
+ cout << "Slop parameter = " << fSlop << endl;
+ cout << "Fiducial volume test flag = " << fvTest << endl;
+ cout << "Fiducial volume excl. width = " << fvDelta << endl;
+ cout << "Decode debug flag = " << fdbg_decoded_cal << endl;
+ cout << "Sparsify debug flag = " << fdbg_sparsified_cal << endl;
+ cout << "Cluster debug flag = " << fdbg_clusters_cal << endl;
+ cout << "Tracking debug flag = " << fdbg_tracks_cal << endl;
BlockThick = new Double_t [fNLayers];
fNBlocks = new Int_t [fNLayers];
@@ -220,6 +230,18 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
}
+ // Fiducial volume limits, based on Plane 1 positions.
+ //
+
+ fvXmin = XPos[0][0] + fvDelta;
+ fvXmax = XPos[0][fNBlocks[0]-1] + BlockThick[0] - fvDelta;
+ fvYmin = YPos[0] + fvDelta;
+ fvYmax = YPos[1] - fvDelta;
+
+ cout << "Fiducial volume limits:" << endl;
+ cout << " Xmin = " << fvXmin << " Xmax = " << fvXmax << endl;
+ cout << " Ymin = " << fvYmin << " Ymax = " << fvYmax << endl;
+
//Calibration related parameters (from hcal.param).
fNtotBlocks=0; //total number of blocks
@@ -430,8 +452,10 @@ Int_t THcShower::DefineVariables( EMode mode )
{ "mult", "Multiplicity of largest cluster", "fMult" },
// { "nblk", "Numbers of blocks in main cluster", "fNblk" },
// { "eblk", "Energies of blocks in main cluster", "fEblk" },
- // { "trx", "track x-position in det plane", "fTRX" },
- // { "try", "track y-position in det plane", "fTRY" },
+ { "trx", "track x-position in det plane (1st track)", "fTRX" },
+ { "try", "track y-position in det plane (1st track)", "fTRY" },
+ { "tre", "track energy in the calorimeter (1st track)", "fTRE" },
+ { "trepr", "track energy in the preshower (1st track)", "fTREpr" },
{ 0 }
};
return DefineVarsFromList( vars, mode );
@@ -486,7 +510,11 @@ void THcShower::Clear(Option_t* opt)
fMult = 0;
fE = 0.;
fEpr = 0.;
- fX = -75.; //out of acceptance
+ fX = -75.; //out of acceptance
+ fTRX = -75.; //out of acceptance
+ fTRY = -40.; //out of acceptance
+ fTRE = -1.;
+ fTREpr = -1.;
}
//_____________________________________________________________________________
@@ -502,7 +530,9 @@ Int_t THcShower::Decode( const THaEvData& evdata )
Int_t nexthit = 0;
for(Int_t ip=0;ip<fNLayers;ip++) {
nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
- //cout << "nexthit = " << nexthit << endl;
+ if (fdbg_decoded_cal) {
+ cout << "THcShower::Decode: nexthit = " << nexthit << endl;
+ }
}
fAnalyzePedestals = 1; // Analyze pedestals first normal events
return(0);
@@ -521,13 +551,16 @@ Int_t THcShower::Decode( const THaEvData& evdata )
}
// fRawHitList is TClones array of THcRawShowerHit objects
- // cout << "THcShower::Decode: Shower raw hit list:" << endl;
- // for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
- // THcRawShowerHit* hit = (THcRawShowerHit *) fRawHitList->At(ihit);
- // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
- // << hit->fADC_pos << " " << hit->fADC_neg << " " << endl;
+
+ // if (fdbg_decoded_cal) {
+ // cout << "THcShower::Decode: Shower raw hit list:" << endl;
+ // for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
+ // THcRawShowerHit* hit = (THcRawShowerHit *) fRawHitList->At(ihit);
+ // cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
+ // << hit->fADC_pos << " " << hit->fADC_neg << " " << endl;
+ // }
+ // cout << endl;
// }
- // cout << endl;
return nhits;
}
@@ -565,7 +598,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
// Clustering of hits.
//
- THcShowerHitList HitList; //list of unclusterd hits
+ THcShowerHitList HitList; //list of unclustered hits
for(Int_t j=0; j < fNLayers; j++) {
@@ -672,30 +705,54 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
fX = (*MaxCluster).clX();
}
- if (fdbg_clusters_cal)
- cout << fEpr << " " << fE << " " << fX << " PrSh" << endl;
+ if (fdbg_clusters_cal)
+ cout << "Max.cluster: fEpr = " << fEpr << " fE = " << fE << " fX = " << fX
+ << endl;
// Track-to-cluster matching.
//
Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
- cout << "Number of reconstructed tracks = " << Ntracks << endl;
+ if (fdbg_tracks_cal) {
+ cout << endl;
+ cout << "Number of reconstructed tracks = " << Ntracks << endl;
+ }
for (Int_t i=0; i<Ntracks; i++) {
THaTrack* theTrack = static_cast<THaTrack*>( tracks[i] );
- // cout << " Track " << i << ": "
- // << " X = " << theTrack->GetX()
- // << " Y = " << theTrack->GetY()
- // << " Theta = " << theTrack->GetTheta()
- // << " Phi = " << theTrack->GetPhi()
- // << endl;
+ if (fdbg_tracks_cal) {
+ cout << " Track " << i << ": "
+ << " X = " << theTrack->GetX()
+ << " Y = " << theTrack->GetY()
+ << " Theta = " << theTrack->GetTheta()
+ << " Phi = " << theTrack->GetPhi()
+ << endl;
+ }
+
+ Double_t Xtr;
+ Double_t Ytr;
+
+ Int_t mclust = MatchCluster(theTrack, ClusterList, Xtr, Ytr);
+
+ if (i==0) {
+ fTRX = Xtr;
+ fTRY = Ytr;
+ if (mclust >= 0) {
+ THcShowerCluster* cluster = (*ClusterList).ListedCluster(mclust);
+ fTRE = (*cluster).clE();
+ fTREpr = (*cluster).clEpr();
+ }
+ }
- MatchCluster(theTrack);
}
+ if (fdbg_tracks_cal)
+ cout << "1st track cluster: fTRX = " << fTRX << " fTRY = " << fTRY
+ << endl;
+
if (fdbg_clusters_cal)
cout << "THcShower::CoarseProcess return ---------------------------" <<endl;
@@ -704,32 +761,102 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
//-----------------------------------------------------------------------------
-void THcShower::MatchCluster(THaTrack* Track)
+Int_t THcShower::MatchCluster(THaTrack* Track,
+ THcShowerClusterList* ClusterList,
+ Double_t& XTrFront, Double_t& YTrFront)
{
// Match a cluster to a given track.
- cout << "Track at DC:"
- << " X = " << Track->GetX()
- << " Y = " << Track->GetY()
- << " Theta = " << Track->GetTheta()
- << " Phi = " << Track->GetPhi()
- << endl;
+ if (fdbg_tracks_cal) {
+ cout << "THcShower::MatchCluster: Track at DC:"
+ << " X = " << Track->GetX()
+ << " Y = " << Track->GetY()
+ << " Theta = " << Track->GetTheta()
+ << " Phi = " << Track->GetPhi()
+ << endl;
+ }
+
+ XTrFront = kBig;
+ YTrFront = kBig;
+ Double_t pathl = kBig;
+
+ // Track interception with face of the calorimeter. The coordinates are
+ // in the calorimeter's local system.
+
+ fOrigin = fPlanes[0]->GetOrigin();
+
+ CalcTrackIntercept(Track, pathl, XTrFront, YTrFront);
+
+ // Transform coordiantes to the spectrometer's coordinate system.
- Double_t xc = 1.e8;
- Double_t yc = 1.e8;
- Double_t pathl = 1.e8;
+ XTrFront += GetOrigin().X();
+ YTrFront += GetOrigin().Y();
- CalcTrackIntercept(Track, pathl, xc, yc);
+ if (fdbg_tracks_cal)
+ cout << " Track at the front of Calorimeter:"
+ << " X = " << XTrFront
+ << " Y = " << YTrFront
+ << " Pathl = " << pathl
+ << endl;
+
+ Bool_t inFidVol = true;
+
+ if (fvTest) {
+
+ fOrigin = fPlanes[fNLayers-1]->GetOrigin();
+
+ Double_t XTrBack = kBig;
+ Double_t YTrBack = kBig;
+
+ CalcTrackIntercept(Track, pathl, XTrBack, YTrBack);
+
+ XTrBack += GetOrigin().X();
+ YTrBack += GetOrigin().Y();
+
+ if (fdbg_tracks_cal)
+ cout << " Track at the back of Calorimeter:"
+ << " X = " << XTrBack
+ << " Y = " << YTrBack
+ << " Pathl = " << pathl
+ << endl;
+
+ inFidVol = (XTrFront < fvXmax) && (XTrFront > fvXmin) &&
+ (YTrFront < fvYmax) && (YTrFront > fvYmin) &&
+ (XTrBack < fvXmax) && (XTrBack > fvXmin) &&
+ (YTrBack < fvYmax) && (YTrBack > fvYmin);
+
+ if (fdbg_tracks_cal)
+ cout << " Fiducial volume test: inFidVol = " << inFidVol << endl;
+ }
+
+ // Match a cluster to the track.
+
+ Int_t mclust = -1; // The match cluster #, initialize with a bogus value.
+ Double_t deltaX = kBig; // Track to cluster distance
+
+ if (inFidVol) {
+
+ for (Int_t i=0; i<fNclust; i++) {
+
+ THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+
+ Double_t dx = TMath::Abs( (*cluster).clX() - XTrFront );
+
+ if (dx <= (0.5*BlockThick[0] + fSlop)) {
+
+ if (dx <= deltaX) {
+ mclust = i;
+ deltaX = dx;
+ }
+ }
+ }
+ }
- // Double_t dx = 0.;
- // Int_t pad = -1;
- // new ( (*fTrackProj)[0]) THaTrackProj(xc,yc,pathl,dx,pad,this);
+ if (fdbg_tracks_cal)
+ cout << "MatchCluster: mclust= " << mclust << " delatX= " << deltaX
+ << endl;
- cout << "Track at Calorimeter:"
- << " X = " << xc
- << " Y = " << yc
- << " Pathl = " << pathl
- << endl;
+ return mclust;
}
diff --git a/src/THcShower.h b/src/THcShower.h
index 03573ba3998bcf94626a1aece453108d5ed1e0b0..76ae8fd7022d00630bd56c12863d163507036867 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -11,6 +11,7 @@
#include "THaNonTrackingDetector.h"
#include "THcHitList.h"
#include "THcShowerPlane.h"
+#include "THcShowerCluster.h"
class THaScCalib;
@@ -35,6 +36,22 @@ public:
Int_t GetNBlocks(Int_t NLayer) const { return fNBlocks[NLayer];}
+ Double_t GetXPos(Int_t NLayer, Int_t NRaw) const {
+ return XPos[NLayer][NRaw];
+ }
+
+ Double_t GetYPos(Int_t NLayer, Int_t Side) const {
+
+ //Side = 0 for postive (left) side
+ //Side = 1 for negative (right) side
+
+ return YPos[2*NLayer+(1-Side)];
+ }
+
+ Double_t GetZPos(Int_t NLayer) const {return fNLayerZPos[NLayer];}
+
+ Double_t GetBlockThick(Int_t NLayer) {return BlockThick[NLayer];}
+
// friend class THaScCalib; not needed for now.
Int_t GetPedLimit(Int_t NBlock, Int_t NLayer, Int_t Side) {
@@ -90,10 +107,12 @@ public:
Double_t fX; // x-position (cm) of the largest cluster
Double_t fZ; // z-position (cm) of the largest cluster
Int_t fMult; // # of hits in the largest cluster
- // Int_t fNblk; // Number of blocks in main cluster
- // Double_t* fEblk; // Energies of blocks in main cluster
- // Double_t fTRX; // track x-position in det plane"
- // Double_t fTRY; // track y-position in det plane",
+ // Int_t fNblk; // Number of blocks in main cluster
+ // Double_t* fEblk; // Energies of blocks in main cluster
+ Double_t fTRX; // track x-position in det plane (1st track)
+ Double_t fTRY; // track y-position in det plane (1st track)
+ Double_t fTRE; // Energy (MeV) of the cluster associated to track
+ Double_t fTREpr; // Preshower Energy (MeV) of the track's cluster
// Potential Hall C parameters. Mostly here for demonstration
@@ -106,10 +125,21 @@ public:
Double_t** XPos; // [fNLayers] X,Y,Z positions of blocks
Double_t* YPos;
Double_t* ZPos;
- Int_t fNegCols; //number of columns with PMTTs on the negative side only.
- Double_t fSlop; //Track to cluster vertical slop distance.
- Int_t fvTest; //fiducial volume test flag
- Int_t fdbg_clusters_cal; // Shower debug flag
+ Int_t fNegCols; // # of columns with neg. side PMTs only.
+ Double_t fSlop; // Track to cluster vertical slop distance.
+ Int_t fvTest; // fiducial volume test flag for tracking
+ Double_t fvDelta; // Exclusion band width for fiducial volume
+
+ Double_t fvXmin; // Fiducial volume limits
+ Double_t fvXmax;
+ Double_t fvYmin;
+ Double_t fvYmax;
+
+
+ Int_t fdbg_decoded_cal; // Shower debug flags
+ Int_t fdbg_sparsified_cal;
+ Int_t fdbg_clusters_cal;
+ Int_t fdbg_tracks_cal;
THcShowerPlane** fPlanes; // [fNLayers] Shower Plane objects
@@ -122,7 +152,9 @@ public:
void Setup(const char* name, const char* description);
- void MatchCluster(THaTrack*);
+ Int_t MatchCluster(THaTrack*, THcShowerClusterList*, Double_t&, Double_t&);
+
+ friend class THcShowerPlane;
ClassDef(THcShower,0) // Generic class
};
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index d579a85ba14255d58fd50dc36403e6e470b62c4c..509ff3ad9831417f0c15ca69a93cdedace866705 100644
--- a/src/THcShowerPlane.cxx
+++ b/src/THcShowerPlane.cxx
@@ -16,6 +16,8 @@
#include "THcRawShowerHit.h"
#include "TClass.h"
#include "math.h"
+#include "THaTrack.h"
+#include "THaTrackProj.h"
#include <cstring>
#include <cstdio>
@@ -115,6 +117,30 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date )
// cout << "THcShowerPlane::ReadDatabase: fLayerNum=" << fLayerNum
// << " fNelem=" << fNelem << endl;
+ // Origin of the plane.
+
+ Double_t BlockThick = fParent->GetBlockThick(fLayerNum-1);
+
+ Double_t xOrig = (fParent->GetXPos(fLayerNum-1,0) +
+ fParent->GetXPos(fLayerNum-1,fNelem-1))/2 +
+ BlockThick/2;
+
+ Double_t yOrig = (fParent->GetYPos(fLayerNum-1,0) +
+ fParent->GetYPos(fLayerNum-1,1))/2;
+
+ Double_t zOrig = fParent->GetZPos(fLayerNum-1);
+
+ fOrigin.SetXYZ(xOrig, yOrig, zOrig);
+
+ cout << "Origin of Layer " << fLayerNum << ": "
+ << fOrigin.X() << " " << fOrigin.Y() << " " << fOrigin.Z() << endl;
+
+ // Detector axes. Assume no rotation.
+ //
+ // DefineAxes(0.); not for subdetector
+
+ //
+
fA_Pos = new Double_t[fNelem];
fA_Neg = new Double_t[fNelem];
fA_Pos_p = new Double_t[fNelem];
@@ -214,8 +240,26 @@ Int_t THcShowerPlane::CoarseProcess( TClonesArray& tracks )
// HitCount();
+ /*
+ if (THcShower::fdbg_tracks_cal)
cout << "THcShowerPlane::CoarseProcess called ---------------------" << endl;
+ Int_t Ntracks = tracks.GetLast()+1; // Number of reconstructed tracks
+
+ if (THcShower::fdbg_tracks_cal)
+ cout << " Number of reconstructed tracks = " << Ntracks << endl;
+
+ for (Int_t i=0; i<Ntracks; i++) {
+
+ THaTrack* theTrack = static_cast<THaTrack*>( tracks[i] );
+
+ Double_t pathl;
+ Double_t xtrk;
+ Double_t ytrk;
+ CalcTrackIntercept(theTrack, pathl, xtrk, ytrk);
+ }
+ */
+
return 0;
}
@@ -231,6 +275,12 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
// Assumes that the hit list is sorted by layer, so we stop when the
// plane doesn't agree and return the index for the next hit.
+ THcShower* fParent;
+ fParent = (THcShower*) GetParent();
+
+ if (fParent->fdbg_decoded_cal)
+ cout << "THcShowerPlane::ProcessHits called ----" << endl;
+
Int_t nPosADCHits=0;
Int_t nNegADCHits=0;
fPosADCHits->Clear();
@@ -251,12 +301,19 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
Int_t nrawhits = rawhits->GetLast()+1;
Int_t ihit = nexthit;
- //cout << "nrawhits = " << nrawhits << endl;
- //cout << "nexthit = " << nexthit << endl;
+
+ if (fParent->fdbg_decoded_cal)
+ {
+ cout << " nrawhits = " << nrawhits << endl;
+ cout << " nexthit = " << nexthit << endl;
+ }
+
while(ihit < nrawhits) {
THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
- //cout << "fplane = " << hit->fPlane << " Num = " << fLayerNum << endl;
+ if (fParent->fdbg_decoded_cal)
+ cout << " fplane = " << hit->fPlane << " Num = " << fLayerNum << endl;
+
if(hit->fPlane > fLayerNum) {
break;
}
@@ -268,13 +325,11 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
fA_Pos_p[hit->fCounter-1] = hit->fADC_pos - fPosPed[hit->fCounter -1];
fA_Neg_p[hit->fCounter-1] = hit->fADC_neg - fNegPed[hit->fCounter -1];
- THcShower* fParent;
- fParent = (THcShower*) GetParent();
-
Double_t thresh_pos = fPosThresh[hit->fCounter -1];
if(hit->fADC_pos > thresh_pos) {
- THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+ THcSignalHit *sighit =
+ (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
sighit->Set(hit->fCounter, hit->fADC_pos);
fEpos[hit->fCounter-1] += fA_Pos_p[hit->fCounter-1]*
@@ -284,7 +339,8 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
Double_t thresh_neg = fNegThresh[hit->fCounter -1];
if(hit->fADC_neg > thresh_neg) {
- THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
+ THcSignalHit *sighit =
+ (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
sighit->Set(hit->fCounter, hit->fADC_neg);
fEneg[hit->fCounter-1] += fA_Neg_p[hit->fCounter-1]*
@@ -296,6 +352,10 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
ihit++;
}
+
+ if (fParent->fdbg_decoded_cal)
+ cout << "THcShowerPlane::ProcessHits return ----" << endl;
+
return(ihit);
}
@@ -306,7 +366,10 @@ Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
// arrays for calculating pedestals.
Int_t nrawhits = rawhits->GetLast()+1;
- // cout << "THcScintillatorPlane::AcculatePedestals " << fLayerNum << " " << nexthit << "/" << nrawhits << endl;
+
+ // cout << "THcScintillatorPlane::AcculatePedestals " << fLayerNum << " "
+ // << nexthit << "/" << nrawhits << endl;
+
Int_t ihit = nexthit;
while(ihit < nrawhits) {
THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
@@ -362,8 +425,8 @@ void THcShowerPlane::CalculatePedestals( )
- fNegPed[i]*fNegPed[i]);
fNegThresh[i] = fNegPed[i] + TMath::Min(50., TMath::Max(10., 3.*fNegSig[i]));
- cout << "Ped&Thr: " << fPosPed[i] << " " << fPosThresh[i] << " " <<
- fNegPed[i] << " " << fNegThresh[i] << " " << i+1 << endl;
+ // cout << "Ped&Thr: " << fPosPed[i] << " " << fPosThresh[i] << " " <<
+ // fNegPed[i] << " " << fNegThresh[i] << " " << i+1 << endl;
}
// cout << " " << endl;
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index 9f2467c0d669937c685729fa8cbbe0df8bf5843e..6c0aa5124fc13cfa0cc8c62fc38c44bcec20179d 100644
--- a/src/THcShowerPlane.h
+++ b/src/THcShowerPlane.h
@@ -46,6 +46,10 @@ class THcShowerPlane : public THaSubDetector {
TClonesArray* fParentHitList;
+ TVector3 GetOrigin() {
+ return fOrigin;
+ }
+
Double_t GetEplane() {
return fEplane;
};