diff --git a/podd b/podd index 5dbaea975ff82d567c1ec4ef866634dc97bb55b3..6dc6a0d2995aa24855eaa4f07f7767ed7486be6f 160000 --- a/podd +++ b/podd @@ -1 +1 @@ -Subproject commit 5dbaea975ff82d567c1ec4ef866634dc97bb55b3 +Subproject commit 6dc6a0d2995aa24855eaa4f07f7767ed7486be6f diff --git a/src/THcShower.cxx b/src/THcShower.cxx index e53a9e3cb1bd226927165929d33c3cf486aadda9..14da4dc859e165676ac3fb3de3a8aaaf34904a30 100644 --- a/src/THcShower.cxx +++ b/src/THcShower.cxx @@ -39,7 +39,6 @@ THcShower::THcShower( const char* name, const char* description, THaNonTrackingDetector(name,description,apparatus) { // Constructor -// fTrackProj = new TClonesArray( "THaTrackProj", 5 ); fNLayers = 0; // No layers until we make them } @@ -50,8 +49,10 @@ THcShower::THcShower( ) : // Constructor } +//_____________________________________________________________________________ void THcShower::Setup(const char* name, const char* description) { + cout << "THcShower::Setup called " << GetName() << endl; char prefix[2]; @@ -65,8 +66,6 @@ void THcShower::Setup(const char* name, const char* description) {0} }; - cout << "THcShower::Setup called " << GetName() << endl; - gHcParms->LoadParmValues((DBRequest*)&list,prefix); cout << layernamelist << endl; cout << "Shower Counter: " << fNLayers << " layers" << endl; @@ -95,6 +94,7 @@ void THcShower::Setup(const char* name, const char* description) strcat(desc, fLayerNames[i]); fPlanes[i] = new THcShowerPlane(fLayerNames[i], desc, i+1, this); + cout << "Created Shower Plane " << fLayerNames[i] << ", " << desc << endl; } @@ -276,15 +276,15 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) //Read in parameters from hcal.param Double_t hcal_pos_cal_const[fNtotBlocks]; - // Double_t hcal_pos_gain_ini[fNtotBlocks]; - // Double_t hcal_pos_gain_cur[fNtotBlocks]; - // Int_t hcal_pos_ped_limit[fNtotBlocks]; + // Double_t hcal_pos_gain_ini[fNtotBlocks]; not used + // Double_t hcal_pos_gain_cur[fNtotBlocks]; not used + // Int_t hcal_pos_ped_limit[fNtotBlocks]; not used Double_t hcal_pos_gain_cor[fNtotBlocks]; Double_t hcal_neg_cal_const[fNtotBlocks]; - // Double_t hcal_neg_gain_ini[fNtotBlocks]; - // Double_t hcal_neg_gain_cur[fNtotBlocks]; - // Int_t hcal_neg_ped_limit[fNtotBlocks]; + // Double_t hcal_neg_gain_ini[fNtotBlocks]; not used + // Double_t hcal_neg_gain_cur[fNtotBlocks]; not used + // Int_t hcal_neg_ped_limit[fNtotBlocks]; not used Double_t hcal_neg_gain_cor[fNtotBlocks]; DBRequest list[]={ @@ -303,8 +303,6 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) }; gHcParms->LoadParmValues((DBRequest*)&list, prefix); - //+++ - cout << "hcal_pos_cal_const:" << endl; for (Int_t j=0; j<fNLayers; j++) { for (Int_t i=0; i<fNBlocks[j]; i++) { @@ -313,22 +311,6 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) cout << endl; }; - // cout << "hcal_pos_gain_ini:" << endl; - // for (Int_t j=0; j<fNLayers; j++) { - // for (Int_t i=0; i<fNBlocks[j]; i++) { - // cout << hcal_pos_gain_ini[j*fNBlocks[j]+i] << " "; - // }; - // cout << endl; - // }; - - // cout << "hcal_pos_gain_cur:" << endl; - // for (Int_t j=0; j<fNLayers; j++) { - // for (Int_t i=0; i<fNBlocks[j]; i++) { - // cout << hcal_pos_gain_cur[j*fNBlocks[j]+i] << " "; - // }; - // cout << endl; - // }; - cout << "fShPosPedLimit:" << endl; for (Int_t j=0; j<fNLayers; j++) { for (Int_t i=0; i<fNBlocks[j]; i++) { @@ -345,8 +327,6 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) cout << endl; }; - //--- - cout << "hcal_neg_cal_const:" << endl; for (Int_t j=0; j<fNLayers; j++) { for (Int_t i=0; i<fNBlocks[j]; i++) { @@ -355,22 +335,6 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) cout << endl; }; - // cout << "hcal_neg_gain_ini:" << endl; - // for (Int_t j=0; j<fNLayers; j++) { - // for (Int_t i=0; i<fNBlocks[j]; i++) { - // cout << hcal_neg_gain_ini[j*fNBlocks[j]+i] << " "; - // }; - // // cout << endl; - // }; - - // cout << "hcal_neg_gain_cur:" << endl; - // for (Int_t j=0; j<fNLayers; j++) { - // for (Int_t i=0; i<fNBlocks[j]; i++) { - // cout << hcal_neg_gain_cur[j*fNBlocks[j]+i] << " "; - // }; - // cout << endl; - // }; - cout << "fShNegPedLimit:" << endl; for (Int_t j=0; j<fNLayers; j++) { for (Int_t i=0; i<fNBlocks[j]; i++) { @@ -387,7 +351,7 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) cout << endl; }; - //Calibration constants in GeV per ADC channel. + // Calibration constants (GeV / ADC channel). for (Int_t i=0; i<fNtotBlocks; i++) { fPosGain[i] = hcal_pos_cal_const[i] * hcal_pos_gain_cor[i]; @@ -416,19 +380,13 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) fTREpl_pos_cor = new Double_t [fNLayers]; fTREpl_neg_cor = new Double_t [fNLayers]; - // Origin of the calorimeter, at tha middle of the face of the detector, - // or at the middle of the front plane of the 1-st layer. + // Origin of the calorimeter, at the centre of the face of the detector, + // or at the centre of the front of the 1-st layer. // Double_t xOrig = (XPos[0][0] + XPos[0][fNBlocks[0]-1])/2 + BlockThick[0]/2; Double_t yOrig = (YPos[0] + YPos[1])/2; Double_t zOrig = fNLayerZPos[0]; - // cout << "Origin of the Calorimeter to be set:" << endl; - // cout << " Xorig = " << xOrig << endl; - // cout << " Yorig = " << yOrig << endl; - // cout << " Zorig = " << zOrig << endl; - // << endl; - fOrigin.SetXYZ(xOrig, yOrig, zOrig); cout << "Origin of the Calorimeter:" << endl; @@ -447,7 +405,6 @@ Int_t THcShower::ReadDatabase( const TDatime& date ) } - //_____________________________________________________________________________ Int_t THcShower::DefineVariables( EMode mode ) { @@ -461,46 +418,31 @@ Int_t THcShower::DefineVariables( EMode mode ) // Register variables in global list RVarDef vars[] = { - { "nhits", "Number of hits", "fNhits" }, + { "nhits", "Number of hits", "fNhits" }, // { "a", "Raw ADC amplitudes", "fA" }, // { "a_p", "Ped-subtracted ADC amplitudes", "fA_p" }, // { "a_c", "Calibrated ADC amplitudes", "fA_c" }, // { "asum_p", "Sum of ped-subtracted ADCs", "fAsum_p" }, // { "asum_c", "Sum of calibrated ADCs", "fAsum_c" }, - { "nclust", "Number of clusters", "fNclust" }, - { "emax", "Energy of largest cluster", "fE" }, - { "eprmax", "Preshower Energy of largest cluster", "fEpr" }, - { "xmax", "x-position (cm) of largest cluster", "fX" }, + { "nclust", "Number of clusters", "fNclust" }, + { "emax", "Energy of largest cluster", "fE" }, + { "eprmax", "Preshower Energy of largest cluster", "fEpr" }, + { "xmax", "x-position (cm) of largest cluster", "fX" }, // { "z", "z-position (cm) of largest cluster", "fZ" }, - { "mult", "Multiplicity of largest cluster", "fMult" }, + { "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 (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" }, - { "trecor", "Y-corrected track Edep (1st track)", "fTRE_cor" }, - { "treprcor", "Y-corrected track Epr (1st track)", "fTREpr_cor" }, - { "treplcor", "Y-corrected track Edep for planes", "fTREpl_cor" }, + { "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" }, + { "trecor", "Y-corrected track Edep (1st track)", "fTRE_cor" }, + { "treprcor", "Y-corrected track Epr (1st track)", "fTREpr_cor" }, + { "treplcor", "Y-corrected track Edep for planes", "fTREpl_cor" }, { 0 } }; return DefineVarsFromList( vars, mode ); - /* - for (Int_t i=0; i<fNLayers; i++) { - RVarDef vars[] = { - { Form("treplcor%d",i+1), - Form("Y-corrected track Edep for plane %d",i+1), - Form("fTREpl_cor[%d]",i) }, - { 0 } - }; - return DefineVarsFromList( vars, mode ); - }; - */ - - //{ "treplposcor","Y-corrected track pos. Edep per plane", "fTREpl_pos_cor"}, - //{ "treplnegcor","Y-corrected track neg. Edep per plane", "fTREpl_neg_cor"}, - return kOK; } @@ -600,6 +542,8 @@ Int_t THcShower::Decode( const THaEvData& evdata ) // fRawHitList is TClones array of THcRawShowerHit objects + // This debug output does not work, needs to be corrected. + // if (fdbg_decoded_cal) { // cout << "THcShower::Decode: Shower raw hit list:" << endl; // for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) { @@ -634,38 +578,36 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks) // reconstructed in THaVDC::CoarseTrack() are used. // // Apply corrections and reconstruct the complete hits. - // - // static const Double_t sqrt2 = TMath::Sqrt(2.); if (fdbg_clusters_cal) cout << "THcShower::CoarseProcess called ---------------------------" <<endl; - // ApplyCorrections(); + // Clustering of hits. + // - // - // Clustering of hits. - // + // Fill list of unclustered hits. + + if (fdbg_clusters_cal) cout << "Filling list of unclustered hits..." << endl; - THcShowerHitList HitList; //list of unclustered hits + THcShowerHitList HitList; for(Int_t j=0; j < fNLayers; j++) { - //cout << "Plane " << j << " Eplane = " << fPlanes[j]->GetEplane() << endl; + if (fdbg_clusters_cal)cout << " Plane " << j << " Eplane = " + << fPlanes[j]->GetEplane() << endl; for (Int_t i=0; i<fNBlocks[j]; i++) { //May be should be done this way. // // Double_t Edep = fPlanes[j]->GetEmean(i); - // if (Edep > 0.) { //hit - // Double_t x = YPos[j][i] + BlockThick[j]/2.; //top + thick/2 - // Double_t z = fNLayerZPos[j] + BlockThick[j]/2.; //front + thick/2 + // if (Edep > 0.) { //hit + // Double_t x = YPos[j][i] + BlockThick[j]/2.; //top + thick/2 + // Double_t z = fNLayerZPos[j] + BlockThick[j]/2.;//front + thick/2 // THcShowerHit* hit = new THcShowerHit(i,j,x,z,Edep); //ENGINE way. // - // if (fPlanes[j]->GetApos(i)> fPlanes[j]->GetPosThr(i) || - // fPlanes[j]->GetAneg(i)> fPlanes[j]->GetNegThr(i)) { //hit if (fPlanes[j]->GetApos(i) - fPlanes[j]->GetPosPed(i) > fPlanes[j]->GetPosThr(i) - fPlanes[j]->GetPosPed(i) || fPlanes[j]->GetAneg(i) - fPlanes[j]->GetNegPed(i) > @@ -679,15 +621,15 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks) HitList.push_back(hit); - //cout << "Hit: Edep = " << Edep << " X = " << x << " Z = " << z << - // " Block " << i << " Layer " << j << endl; + if (fdbg_clusters_cal) + cout << " Hit: Edep = " << Edep << " X = " << x << " Z = " << z + << " Block " << i << " Layer " << j << endl; }; } } //Print out hits before clustering. - // fNhits = HitList.size(); @@ -699,13 +641,14 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks) } } + // Create list of clusters and fill it. + THcShowerClusterList* ClusterList = new THcShowerClusterList; ClusterList->ClusterHits(HitList); - fNclust = (*ClusterList).NbClusters(); + fNclust = (*ClusterList).NbClusters(); //number of clusters //Print out the cluster list. - // if (fdbg_clusters_cal) { @@ -731,9 +674,10 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks) } } + // The following is for testing the code, by comparison with Engine. - //The largest cluster. - // + // The largest cluster by deposited energy, its size and energy + // depositions. if (fNclust != 0) { @@ -753,9 +697,9 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks) } } - fMult = (*MaxCluster).clSize(); - fEpr = (*MaxCluster).clEpr(); - fX = (*MaxCluster).clX(); + fMult = (*MaxCluster).clSize(); // number of hit counters + fEpr = (*MaxCluster).clEpr(); // Preshower energy + fX = (*MaxCluster).clX(); // X coordinate (in vertical direction) } if (fdbg_clusters_cal) @@ -785,27 +729,39 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks) << endl; } + // Track coordinates at front of the calorimeter, initialize out of + // acceptance. Double_t Xtr = -75.; Double_t Ytr = -40.; + // Associate a cluster to the track. + Int_t mclust = MatchCluster(theTrack, ClusterList, Xtr, Ytr); + // Do this for the 1-st track only for now. + // if (itrk==0) { - fTRX = Xtr; + fTRX = Xtr; // track coordinates at the calorimeter fTRY = Ytr; - if (mclust >= 0) { + if (mclust >= 0) { // if there is a matched cluster + //Assign energies (not Y corrected) of the matched cluster to the track. THcShowerCluster* cluster = (*ClusterList).ListedCluster(mclust); fTRE = (*cluster).clE(); fTREpr = (*cluster).clEpr(); - fTRE_cor = 0.; + // Correct track energy depositions for the impact coordinate. + + fTRE_cor = 0.; // Coord. corrected total energy depostion for (Int_t ip=0; ip<fNLayers; ip++) { - Float_t corpos = 1.; + // Coordinate correction factors for positive and negative sides, + // different for single PMT counters in the 1-st two layes and for + // 2 PMT counters in the rear two layers. + Float_t corpos = 1.; Float_t corneg = 1.; if (ip < fNegCols) { corpos = Ycor(Ytr,0); @@ -835,7 +791,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks) } //over planes - fTREpr_cor = fTREpl_cor[0]; + fTREpr_cor = fTREpl_cor[0]; // Coord. corrected Edep in Preshower if (fdbg_tracks_cal) cout << " fTREpr = " << fTREpr @@ -863,7 +819,8 @@ Int_t THcShower::MatchCluster(THaTrack* Track, THcShowerClusterList* ClusterList, Double_t& XTrFront, Double_t& YTrFront) { - // Match a cluster to a given track. + // Match a cluster to a given track. Return the cluster number, + // and track coordinates at the front of calorimeter. if (fdbg_tracks_cal) { cout << "THcShower::MatchCluster: Track at DC:" @@ -897,10 +854,15 @@ Int_t THcShower::MatchCluster(THaTrack* Track, << " Pathl = " << pathl << endl; - Bool_t inFidVol = true; + Bool_t inFidVol = true; // In Fiducial Volume flag + + // Re-evaluate Fid. Volume Flag if fid. volume test is requested if (fvTest) { + // Track coordinates at the back of the detector. + + // Origin at the front of the last layer. fOrigin = fPlanes[fNLayers-1]->GetOrigin(); Double_t XTrBack = kBig; @@ -908,8 +870,8 @@ Int_t THcShower::MatchCluster(THaTrack* Track, CalcTrackIntercept(Track, pathl, XTrBack, YTrBack); - XTrBack += GetOrigin().X(); - YTrBack += GetOrigin().Y(); + XTrBack += GetOrigin().X(); // from local coord. system + YTrBack += GetOrigin().Y(); // to the spectrometer system if (fdbg_tracks_cal) cout << " Track at the back of Calorimeter:" diff --git a/src/THcShower.h b/src/THcShower.h index 0452731cee3a76e65b4c7bf2034b8719b06b0164..4286a8252b1b169efb2b52fd473c49adf22a7b44 100644 --- a/src/THcShower.h +++ b/src/THcShower.h @@ -32,8 +32,9 @@ public: Int_t GetNHits() const { return fNhits; } - Int_t GetNTracks() const { return fTrackProj->GetLast()+1; } - const TClonesArray* GetTrackHits() const { return fTrackProj; } + // These are not used in the calorimeter code. + // Int_t GetNTracks() const { return fTrackProj->GetLast()+1; } + // const TClonesArray* GetTrackHits() const { return fTrackProj; } Int_t GetNBlocks(Int_t NLayer) const { return fNBlocks[NLayer];} @@ -43,8 +44,8 @@ public: Double_t GetYPos(Int_t NLayer, Int_t Side) const { - //Side = 0 for postive (left) side - //Side = 1 for negative (right) side + //Side = 0 for postive (right) side + //Side = 1 for negative (left) side return YPos[2*NLayer+(1-Side)]; } @@ -90,11 +91,10 @@ public: } //Coordinate correction for single PMT modules. - //PMT attached at right (negative) side. + //PMT attached at right (positive) side. Float_t Ycor(Double_t y) { return TMath::Exp(y/fAcor)/(1. + y*y/fBcor); - // return TMath::Exp(-y/fAcor)/(1. + y*y/fBcor); } //Coordinate correction for double PMT modules. @@ -106,7 +106,6 @@ public: return 0.; } Int_t sign = 1 - 2*side; - // Int_t sign = 2*side - 1; return (fCcor + sign*y)/(fCcor + sign*y/fDcor); } @@ -135,9 +134,13 @@ protected: 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 + + // Track quantities, are here for test purposes. Better to move + // to the track class(es) later on. + 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 of the cluster associated to track + Double_t fTRE; // Energy of a cluster associated to the track Double_t fTREpr; // Preshower Energy of the track's cluster Double_t fTRE_cor; // Y-corrected track total energy Double_t fTREpr_cor; // Y-corrected track Preshower energy @@ -193,7 +196,6 @@ protected: Int_t MatchCluster(THaTrack*, THcShowerClusterList*, Double_t&, Double_t&); friend class THcShowerPlane; //to access debug flags. - // friend class THcShowerCluster; //to access fNLayers ClassDef(THcShower,0) // Generic class }; diff --git a/src/THcShowerCluster.h b/src/THcShowerCluster.h index 5705abbf024b602660f1ddeea88d81fec664cea5..23d1e0d642cef1ad290c20bee8a89156af69a535 100644 --- a/src/THcShowerCluster.h +++ b/src/THcShowerCluster.h @@ -42,8 +42,8 @@ class THcShowerCluster : THcShowerHitList { (*(THcShowerHitList::begin()+num))->show(); } - // X coordinate of cluster's center of gravity, - // calculated as a weighted by hit energies average. + // X coordinate of center of gravity of cluster, + // calculated as hit energy weighted average. // Put X out of the calorimeter (-75 cm), if there is no energy deposition // in the cluster. // @@ -59,8 +59,8 @@ class THcShowerCluster : THcShowerHitList { return (Etot != 0. ? x_sum/Etot : -75.); } - // Z coordinate of cluster's center of gravity, - // calculated as a weighted by hit energies average. + // Z coordinate of center of gravity of cluster, + // calculated as a hit energy weighted average. // Put Z out of the calorimeter (0 cm), if there is no energy deposition // in the cluster. // @@ -101,10 +101,10 @@ class THcShowerCluster : THcShowerHitList { return Epr; } - //Cluster energy deposition in plane iplane=0,..,3. - // side=0 -- from positive PMTs only - // side=1 -- from negative PMTs only - // side=2 -- from postive and negative PMTs + //Cluster energy deposition in plane iplane=0,..,3: + // side=0 -- from positive PMTs only; + // side=1 -- from negative PMTs only; + // side=2 -- from postive and negative PMTs. // Double_t clEplane(Int_t iplane, Int_t side) { @@ -191,18 +191,18 @@ class THcShowerClusterList : private THcShClusterList { void ClusterHits(THcShowerHitList HitList) { -// Cluster hits from the HitList. The resultant clusters of hits are saved -//in the ClusterList. +// Collect hits from the HitList into the clusters. The resultant clusters +// of hits are saved in the ClusterList. while (HitList.size() != 0) { THcShowerCluster* cluster = new THcShowerCluster; - (*cluster).grow(*(HitList.end()-1)); //move the last hit from the hit list + (*cluster).grow(*(HitList.end()-1)); //Move the last hit from the hit list HitList.erase(HitList.end()-1); //into the 1st cluster bool clustered = true; - while (clustered) { //while a hit is clustered + while (clustered) { //Proceed while a hit is clustered clustered = false; diff --git a/src/THcShowerHit.h b/src/THcShowerHit.h index df7e98feed91365a96c0ae40327ceda04904f3f1..269ddde61836ce95cef507539c478e465fadbdb6 100644 --- a/src/THcShowerHit.h +++ b/src/THcShowerHit.h @@ -12,7 +12,7 @@ using namespace std; class THcShowerHit { //HMS calorimeter hit class private: - Int_t fCol, fRow; //hit colomn and row + Int_t fCol, fRow; //hit colomn and row Double_t fX, fZ; //hit X (vert.) and Z (along spect.axis) coordinates Double_t fE; //hit mean energy deposition Double_t fEpos; //hit energy deposition from positive PMT diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx index bdca666c4e1d3478dbc6af4551430f6c56c217a9..be062a012139ddaa08d09aafa63d01018ac4bfb5 100644 --- a/src/THcShowerPlane.cxx +++ b/src/THcShowerPlane.cxx @@ -116,7 +116,12 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date ) // cout << "THcShowerPlane::ReadDatabase: fLayerNum=" << fLayerNum // << " fNelem=" << fNelem << endl; - // Origin of the plane. + // Origin of the plane: + // + // X is average of top X coordinates of the top and bottom blocks + // shifted by a half of the block thickness; + // Y is average of left and right edges; + // Z is _front_ position of the plane along the beam. Double_t BlockThick = fParent->GetBlockThick(fLayerNum-1); @@ -136,20 +141,9 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date ) // 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]; - fA_Neg_p = new Double_t[fNelem]; + // DefineAxes(0.); Do not work for subdetector - fEpos = new Double_t[fNelem]; - fEneg = new Double_t[fNelem]; - fEmean= new Double_t[fNelem]; - - // fNelem = *(Int_t *)gHcParms->Find(parname)->GetValuePointer(); +// fNelem = *(Int_t *)gHcParms->Find(parname)->GetValuePointer(); // // parname[plen]='\0'; // strcat(parname,"_spacing"); @@ -164,6 +158,8 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date ) // Create arrays to hold results here + // Pedestal limits per channel. + fPosPedLimit = new Int_t [fNelem]; fNegPedLimit = new Int_t [fNelem]; @@ -184,6 +180,19 @@ Int_t THcShowerPlane::ReadDatabase( const TDatime& date ) InitializePedestals(); + // ADC amplitudes per channel. + + fA_Pos = new Double_t[fNelem]; + fA_Neg = new Double_t[fNelem]; + fA_Pos_p = new Double_t[fNelem]; + fA_Neg_p = new Double_t[fNelem]; + + // Energy depositions per block (not corrected for track coordinate) + + fEpos = new Double_t[fNelem]; + fEneg = new Double_t[fNelem]; + fEmean= new Double_t[fNelem]; + return kOK; } @@ -201,14 +210,14 @@ Int_t THcShowerPlane::DefineVariables( EMode mode ) RVarDef vars[] = { {"posadchits", "List of Positive ADC hits","fPosADCHits.THcSignalHit.GetPaddleNumber()"}, {"negadchits", "List of Negative ADC hits","fNegADCHits.THcSignalHit.GetPaddleNumber()"}, - {"apos", "Raw Positive ADC Amplitudes", "fA_Pos"}, - {"aneg", "Raw Negative ADC Amplitudes", "fA_Neg"}, - {"apos_p", "Ped-subtracted Positive ADC Amplitudes", "fA_Pos_p"}, - {"aneg_p", "Ped-subtracted Negative ADC Amplitudes", "fA_Neg_p"}, - {"epos", "Energy Depositions from Positive Side PMTs", "fEpos"}, - {"eneg", "Energy Depositions from Negative Side PMTs", "fEneg"}, - {"emean", "Mean Energy Depositions", "fEmean"}, - {"eplane", "Energy Deposition per plane", "fEplane"}, + {"apos", "Raw Positive ADC Amplitudes", "fA_Pos"}, + {"aneg", "Raw Negative ADC Amplitudes", "fA_Neg"}, + {"apos_p", "Ped-subtracted Positive ADC Amplitudes", "fA_Pos_p"}, + {"aneg_p", "Ped-subtracted Negative ADC Amplitudes", "fA_Neg_p"}, + {"epos", "Energy Depositions from Positive Side PMTs", "fEpos"}, + {"eneg", "Energy Depositions from Negative Side PMTs", "fEneg"}, + {"emean", "Mean Energy Depositions", "fEmean"}, + {"eplane", "Energy Deposition per plane", "fEplane"}, {"eplane_pos", "Energy Deposition per plane from pos. PMTs","fEplane_pos"}, {"eplane_neg", "Energy Deposition per plane from neg. PMTs","fEplane_neg"}, { 0 } @@ -238,7 +247,10 @@ Int_t THcShowerPlane::Decode( const THaEvData& evdata ) //_____________________________________________________________________________ Int_t THcShowerPlane::CoarseProcess( TClonesArray& tracks ) { - + + // Nothing is done here. See ProcessHits method instead. + // + // HitCount(); /* @@ -270,6 +282,7 @@ Int_t THcShowerPlane::FineProcess( TClonesArray& tracks ) return 0; } +//_____________________________________________________________________________ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) { // Extract the data for this layer from hit list @@ -282,6 +295,8 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) if (fParent->fdbg_decoded_cal) cout << "THcShowerPlane::ProcessHits called ----" << endl; + // Initialize variable. + Int_t nPosADCHits=0; Int_t nNegADCHits=0; fPosADCHits->Clear(); @@ -301,6 +316,9 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) fEplane_pos = 0; fEplane_neg = 0; + // Process raw hits. Get ADC hits for the plane, assign variables for each + // channel. + Int_t nrawhits = rawhits->GetLast()+1; Int_t ihit = nexthit; @@ -317,6 +335,8 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) if (fParent->fdbg_decoded_cal) cout << " fplane = " << hit->fPlane << " Num = " << fLayerNum << endl; + // This is OK as far as the hit list is sorted by layer. + // if(hit->fPlane > fLayerNum) { break; } @@ -328,6 +348,9 @@ 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]; + // Sparsify positive side hits, fill the hit list, compute the + // energy depostion from positive side for the counter. + Double_t thresh_pos = fPosThresh[hit->fCounter -1]; if(hit->fADC_pos > thresh_pos) { @@ -339,6 +362,9 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) fParent->GetGain(hit->fCounter-1,fLayerNum-1,0); } + // Sparsify negative side hits, fill the hit list, compute the + // energy depostion from negative side for the counter. + Double_t thresh_neg = fNegThresh[hit->fCounter -1]; if(hit->fADC_neg > thresh_neg) { @@ -350,7 +376,12 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit) fParent->GetGain(hit->fCounter-1,fLayerNum-1,1); } + // Mean energy in the counter. + fEmean[hit->fCounter-1] += (fEpos[hit->fCounter-1] + fEneg[hit->fCounter-1]); + + // Accumulate energies in the plane. + fEplane += fEmean[hit->fCounter-1]; fEplane_pos += fEpos[hit->fCounter-1]; fEplane_neg += fEneg[hit->fCounter-1]; @@ -377,8 +408,12 @@ Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit) Int_t ihit = nexthit; while(ihit < nrawhits) { + THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit); + //cout << "fPlane = " << hit->fPlane << " Limit = " << fLayerNum << endl; + + // OK for hit list sorted by layer. if(hit->fPlane > fLayerNum) { break; } diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h index 8b4de6f8e83fc478b16d116bf85ff1f1011033f4..9bf2388a3975f2327f2b287a800689911fa6048e 100644 --- a/src/THcShowerPlane.h +++ b/src/THcShowerPlane.h @@ -42,7 +42,7 @@ public: virtual Int_t AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit); virtual void CalculatePedestals( ); - Double_t fSpacing; + // Double_t fSpacing; not used TClonesArray* fParentHitList; @@ -120,8 +120,8 @@ protected: Double_t fEplane_pos; // Energy deposition in the plane from positive PMTs Double_t fEplane_neg; // Energy deposition in the plane from negative PMTs - - TClonesArray* fPosADCHits; // List of positive ADC hits + // These lists are not used actively for now. + TClonesArray* fPosADCHits; // List of positive ADC hits TClonesArray* fNegADCHits; // List of negative ADC hits Int_t fLayerNum; // Layer # 1-4