diff --git a/TCatTpcHitClusterProcessor.cc b/TCatTpcHitClusterProcessor.cc index 90416f487767ca47ec67d1c8e62096d7680a222c..1a7aaae29aeb2c8ad2715317f20d3473a1ad481e 100644 --- a/TCatTpcHitClusterProcessor.cc +++ b/TCatTpcHitClusterProcessor.cc @@ -3,7 +3,7 @@ * @brief make cluster of tpc hits * * @date Created : 2016-04-17 04:46:27 JST - * Last Modified : 2019-07-04 21:59:25 JST (ota) + * Last Modified : 2020-05-20 11:04:41 JST (ota) * @author Shinsuke OTA * * (C) 2016 Shinsuke OTA @@ -23,6 +23,7 @@ #include #include #include +#include "TTpcStaticGasProperty.h" using art::TCatTpcHitClusterProcessor; ClassImp(TCatTpcHitClusterProcessor) @@ -74,7 +75,6 @@ namespace { } TCatTpcHitClusterProcessor::TCatTpcHitClusterProcessor() -: fInputName(""), fInput(NULL), fOutputName(""), fOutput(NULL), fReadoutPadArrayName(""), fReadoutPadArray(NULL) { DoubleVec_t defClusterRadius; @@ -84,26 +84,22 @@ TCatTpcHitClusterProcessor::TCatTpcHitClusterProcessor() defPulseLength[0] = 0; defPulseLength[1] = 100; - RegisterInputCollection("InputCollection","input pulse collection (for a specific region)",fInputName,TString("atpulse"), - &fInput,TClonesArray::Class_Name(),TCatPulseShape::Class_Name()); - RegisterOutputCollection("OutputCollection","output pulse collection as a candidate of recoil track",fOutputName,TString("atcluster"), - &fOutput,TClonesArray::Class_Name(),TCatTpcTrack::Class_Name()); - RegisterInputInfo("ReadoutPad","Name of readout pad array",fReadoutPadArrayName,TString("readoutPad"), - &fReadoutPadArray,TCatReadoutPadArray::Class_Name()); - RegisterProcessorParameter("ClusterRadius","radius of cluster",fClusterRadius,defClusterRadius); - RegisterProcessorParameter("MinClusterSize","required minimum cluster size ",fMinClusterSize,5.); - RegisterProcessorParameter("IsSingle","select maximum size cluster if 1 (default: 0)",fIsSingle,kTRUE); - RegisterProcessorParameter("MinHeight","minimum pulse height",fMinHeight,0.); - RegisterProcessorParameter("PulseLengthRange","range of pulse length [start,end]. ([0,100])",fPulseLengthRange,defPulseLength); - RegisterProcessorParameter("SelectIf","select all if 0, inside if 1 and outside if 2. The region is indicated by SelectedRegion", - fSelectIf,0); - RegisterProcessorParameter("SelectedRegion","selected region defined by box [x1,y1,x2,y2]",fSelectedRegion,DoubleVec_t(4,0.)); - RegisterProcessorParameter("ReadoutPadType","select readout pad type 2 if 2 (default: 2)",fReadoutPadType,2); - RegisterProcessorParameter("DoTrackUpdate","flag for updating track",fDoTrackUpdate,kFALSE); + + + Register(fInput("InputCollection","input pulse collection (for a specific region)",TString("atpulse"))); + Register(fOutput("OutputCollection","output pulse collection as a candidate of recoil track",TString("atcluster"))); + Register(fReadoutPadArray("ReadoutPad","Name of readout pad array",TString("readoutPad"))); + Register(fTpcProperty("TpcProperty","Name of tpc gas property",TString("tpcGasProperty"))); + Register(fClusterRadius("ClusterRadius","radius of cluster",defClusterRadius)); + Register(fMinClusterSize("MinClusterSize","required minimum cluster size ",5.)); + Register(fIsSingle("IsSingle","select maximum size cluster if 1 (default: 0)",kTRUE)); + Register(fMinHeight("MinHeight","minimum pulse height",0.)); + Register(fPulseLengthRange("PulseLengthRange","range of pulse length [start,end]. ([0,100])",defPulseLength)); + Register(fSelectIf("SelectIf","select all if 0, inside if 1 and outside if 2. The region is indicated by SelectedRegion",0)); + Register(fSelectedRegion("SelectedRegion","selected region defined by box [x1,y1,x2,y2]",DoubleVec_t(4,0.))); + Register(fDoTrackUpdate("DoTrackUpdate","flag for updating track",kFALSE)); RegisterHelper(&fTrackUpdateHelper); - fTrackUpdateHelper.SetUpdateTracks(&fOutput); - RegisterProcessorParameter("ChargeWeightPower","order of power of weighted charge",fChargeWeightPower,1.); - RegisterProcessorParameter("DriftVelocity", "drift velocity (cm/us)",fDriftVelocity,1.); + fTrackUpdateHelper.SetUpdateTracks((TClonesArray**)fOutput.PtrRef()); } TCatTpcHitClusterProcessor::~TCatTpcHitClusterProcessor() @@ -127,36 +123,35 @@ void TCatTpcHitClusterProcessor::Init(TEventCollection *col) { //====================================================================== // check the type of readout pad aray - Int_t id = (*fReadoutPadArray)->GetUniqueID() ; -Info("Init","Readout Pad Type = %d",id); + Int_t id = fReadoutPadArray->GetUniqueID() ; + Info("Init","Readout Pad Type = %d",id); //====================================================================== // cluster radii for XZ and Y, respectively - if (fClusterRadius.size() != 2) { - SetStateError(Form("ClusterRadius has %d components while ClusterRadius should have two values for ZX plane and Y direction, respectively\n",(Int_t)fClusterRadius.size())); + if (((DoubleVec_t)fClusterRadius).size() != 2) { + SetStateError(Form("ClusterRadius has %d components while ClusterRadius should have two values for ZX plane and Y direction, respectively\n",(Int_t)((DoubleVec_t)fClusterRadius).size())); return; } //====================================================================== // cluster radii for XZ and Y, respectively + DoubleVec_t &selectedRegion = fSelectedRegion.Value(); if (fSelectIf) { - if (fSelectedRegion.size() != 4) { - SetStateError(Form("SelectedRegion should be defined by box [x1,y1,x2,y2] (size = %d)\n",(Int_t)fSelectedRegion.size())); + if (selectedRegion.size() != 4) { + SetStateError(Form("SelectedRegion should be defined by box [x1,y1,x2,y2] (size = %d)\n",(Int_t)selectedRegion.size())); } - if (fSelectedRegion[0] > fSelectedRegion[2]) { - Double_t tmp = fSelectedRegion[0] ; - fSelectedRegion[0] = fSelectedRegion[2]; - fSelectedRegion[2] = tmp; + if (selectedRegion[0] > selectedRegion[2]) { + Double_t tmp = selectedRegion[0] ; + selectedRegion[0] = selectedRegion[2]; + selectedRegion[2] = tmp; } - if (fSelectedRegion[1] > fSelectedRegion[3]) { - Double_t tmp = fSelectedRegion[1] ; - fSelectedRegion[1] = fSelectedRegion[3]; - fSelectedRegion[3] = tmp; + if (selectedRegion[1] > selectedRegion[3]) { + Double_t tmp = selectedRegion[1] ; + selectedRegion[1] = selectedRegion[3]; + selectedRegion[3] = tmp; } } - - fDriftVelocity *= (TArtSystemOfUnit::cm / TArtSystemOfUnit::us); } void TCatTpcHitClusterProcessor::Process() @@ -167,48 +162,43 @@ void TCatTpcHitClusterProcessor::Process() // initialize local variables pulseVec_t input; clusterVec_t cluster; + const DoubleVec_t& radius = fClusterRadius.Value(); + const DoubleVec_t& reg = fSelectedRegion.Value(); - Int_t nHits = (*fInput)->GetEntriesFast(); + Int_t nHits = fInput->GetEntriesFast(); + std::vector use(nHits,true); TCatPulseShape::SetSortType(TCatPulseShape::kID); TCatPulseShape::SetSortOrder(TCatPulseShape::kASC); - (*fInput)->Sort(); - (*fInput)->Compress(); + fInput->Sort(); + fInput->Compress(); - // ignore beam particle and bad pads + // ignore bad pads, invalid pulse, outside/inside given region for (Int_t iHit = 0; iHit != nHits; iHit++) { - TCatPulseShape *pulse = (TCatPulseShape*) (*fInput)->UncheckedAt(iHit); - TCatReadoutPad *pad = static_cast((*fReadoutPadArray)->UncheckedAt(pulse->GetID())); - if (pad->IsBad()) continue; - - Bool_t doUse = kTRUE; - if (fSelectIf) { - // judge if pad is inside region - Double_t padx = pulse->GetX(); - Double_t pady = pulse->GetZ(); - Double_t isInside = kTRUE; - if (padx < fSelectedRegion[0] || fSelectedRegion[2] < padx || - pady < fSelectedRegion[1] || fSelectedRegion[3] < pady) { - // outside region - isInside = kFALSE; - } - switch (fSelectIf) { - case 1: - if (!isInside) doUse = kFALSE; - break; - case 2: - if (isInside) doUse = kFALSE; - break; - default: - break; - } - } - - if (!doUse) continue; - input.push_back(pulse); - + TCatPulseShape *pulse = (TCatPulseShape*) fInput->UncheckedAt(iHit); + TCatReadoutPad *pad = static_cast(fReadoutPadArray->UncheckedAt(pulse->GetID())); + if (!pulse->IsValid() || pad->IsBad()) { + use[iHit] = false; + continue; + } + if (!fSelectIf) { + // assume use is initialized to be true + continue; + } + const Double_t padx = pulse->GetX(); + const Double_t pady = pulse->GetZ(); + const bool isInside = reg[0] < padx && padx < reg[2] && reg[1] < pady && pady < reg[3] ; + use[iHit] = ((fSelectIf == 1) && isInside) || ((fSelectIf == 2) && !isInside); } + for (Int_t iHit = 0; iHit != nHits; iHit++) { + if (!use[iHit]) continue; + input.push_back((TCatPulseShape*) fInput->UncheckedAt(iHit)); + } + +// printf("input size = %d\n",input.size()); + + for (pulseVec_t::iterator it1 = input.begin(), itend = input.end(); it1 != itend; ++it1) { Bool_t used = kFALSE; for (clusterVec_t::iterator itCluster = cluster.begin(), itClusterEnd = cluster.end(); @@ -218,12 +208,14 @@ void TCatTpcHitClusterProcessor::Process() itHitInCls != itEndHitInCls; ++itHitInCls) { - TVector3 p1 = (*it1)->GetPos(); - TVector3 p2 = (*itHitInCls)->GetPos(); - Double_t y1 = p1.Y(); p1.SetY(0); - Double_t y2 = p2.Y(); p2.SetY(0); - - if ((p1-p2).Mag() < fClusterRadius[0] && TMath::Abs(y1-y2) < fClusterRadius[1]) { + const TVector3& p1 = (*it1)->GetPos(); + const TVector3& p2 = (*itHitInCls)->GetPos(); + Double_t y1 = p1.Y(); + Double_t y2 = p2.Y(); + double dx = p2.X() - p1.X(); + double dz = p2.Z() - p1.Z(); + double mag = TMath::Sqrt(dx * dx + dz * dz); + if (mag < radius[0] && TMath::Abs(y1-y2) < radius[1]) { (*itCluster).push_back(*it1); used = kTRUE; break; @@ -238,7 +230,7 @@ void TCatTpcHitClusterProcessor::Process() } // find the neighbor clusters - concatClusters(cluster,fClusterRadius[0],fClusterRadius[1]); + concatClusters(cluster,radius[0],radius[1]); clusterList_t candidates; for (clusterVec_t::iterator itCluster = cluster.begin(), itClusterEnd = cluster.end(); @@ -254,11 +246,9 @@ void TCatTpcHitClusterProcessor::Process() for (clusterList_t::iterator itc = candidates.begin(), itcEnd = candidates.end(); itc != itcEnd; ++itc,++clusterID) { TCatTpcTrack *track = static_cast(fOutput->ConstructedAt(fOutput->GetEntriesFast())); track->Init(); - track->SetChargeWeight(fChargeWeightPower); - track->EstimateFromHits(*itc); track->SetID(clusterID); - track->SetTimestamp(track->Y() / fDriftVelocity); + track->SetTimestamp(track->Y() / fTpcProperty->GetDriftVelocity()); if (fIsSingle) break; } diff --git a/TCatTpcHitClusterProcessor.h b/TCatTpcHitClusterProcessor.h index 5a4b7d9d6071c168774e4fc174f2aeacc2dc7f93..785357d2297ba9798a6a251e5962310420537776 100644 --- a/TCatTpcHitClusterProcessor.h +++ b/TCatTpcHitClusterProcessor.h @@ -3,7 +3,7 @@ * @brief make cluster of tpc hits * * @date Created : 2016-04-17 04:45:26 JST - * Last Modified : 2018-07-18 13:27:58 JST (ota) + * Last Modified : 2020-05-19 16:14:13 JST (ota) * @author Shinsuke OTA * * (C) 2016 Shinsuke OTA @@ -18,6 +18,9 @@ namespace art { class TCatTpcHitClusterProcessor; class TCatReadoutPadArray; + class TCatPulseShape; + class TCatTpcTrack; + class TTpcStaticGasProperty; } class TVector3; @@ -36,40 +39,21 @@ public: virtual void Init(TEventCollection *col); protected: + InputData fInput; //! input array + OutputData fOutput; + InputInfo fReadoutPadArray; + InputInfo fTpcProperty; + Parameter fClusterRadius; + Parameter fIsSingle; + Parameter fMinClusterSize; + Parameter fMinHeight; + Parameter fPulseLengthRange; + Parameter fSelectIf; + Parameter fSelectedRegion; + Parameter fDoTrackUpdate; - TString fInputName; // name of input collection - TClonesArray **fInput; //! input array - TString fOutputName; // name of output collection - TClonesArray *fOutput; //! output array - TString fReadoutPadArrayName; // name of readout pad array - TCatReadoutPadArray **fReadoutPadArray; //! readout pad array - - TString fNameInertia; // name of output moment of inertia - Double_t *fInertia; //! - TString fNamePrincipalAxis; // name of output principal axis - TVector3 *fPrincipalAxis; //! - TString fNameDeviation; // deviation - TVector3 *fDeviation; //! - - - DoubleVec_t fClusterRadius; // size of cluster (mm) - - Double_t fMinClusterSize; // required minimum cluster size - Bool_t fIsSingle; // select maximum size cluster if 1 (default:0) - Int_t fReadoutPadType; // select readout pad type 2 if 2 (default:2) - - Double_t fMinHeight; - DoubleVec_t fPulseLengthRange; - Int_t fSelectIf; // case for usage of selected region: 0 never, 1 accept, 2 reject - DoubleVec_t fSelectedRegion; // rectangular region for selection - - Bool_t fDoTrackUpdate; TCatTpcTrackUpdateHelper fTrackUpdateHelper; - Double_t fChargeWeightPower; //charge weight power - - Double_t fDriftVelocity; // drift velocity (cm/us) - private: diff --git a/TCatTpcTrack.cc b/TCatTpcTrack.cc index dbb1d3c240fbaa9351efedd6fbd1b1605cb9b986..7902222b6d6a0d7b77e5156e9ea79d14ee5ebd03 100644 --- a/TCatTpcTrack.cc +++ b/TCatTpcTrack.cc @@ -3,7 +3,7 @@ * @brief tpc track * * @date Created : 2016-12-09 09:50:14 JST - * Last Modified : 2020-04-03 16:55:13 JST (ota) + * Last Modified : 2020-05-20 09:55:27 JST (ota) * @author Shinsuke OTA * * (C) 2016 Shinsuke OTA @@ -168,37 +168,46 @@ void TCatTpcTrack::EstimateFromHits(const std::vector& clu // set values fInertia.SetXYZ(eigenval[2],eigenval[1],eigenval[0]); fPrincipalAxis.SetXYZ(eigenvec[0][2],eigenvec[1][2],eigenvec[2][2]); + TVector3 secondAxis(eigenvec[0][1],eigenvec[1][1],eigenvec[2][1]); + TVector3 thirdAxis(eigenvec[0][0],eigenvec[1][0],eigenvec[2][0]); // set value + // local variable for projected length search + TVector3 minPoint[3]; + TVector3 masPoint[3]; + std::vector mindist(3,TMath::Limits::Max()); + std::vector maxdist(3,TMath::Limits::Min()); TVector3 centroid(mx,my,mz); - TVector3 unit = fPrincipalAxis.Unit(); - Double_t min = TMath::Limits::Max(); - Double_t max = TMath::Limits::Min(); - std::vector::const_iterator itStartPoint = cluster.begin(); - std::vector::const_iterator itEndPoint = cluster.begin(); - - - for (std::vector::const_iterator it = cluster.begin(), itend = cluster.end(); it < itend; ++it) { - // estimate start and end point - TCatPulseShape *pulse = *it; + TVector3 unit[3] = {fPrincipalAxis.Unit(), secondAxis.Unit(), thirdAxis.Unit()}; + for (int i = 0; i < 3; ++i) { + Double_t min = TMath::Limits::Max(); + Double_t max = TMath::Limits::Min(); + std::vector::const_iterator itStartPoint = cluster.begin(); + std::vector::const_iterator itEndPoint = cluster.begin(); + for (std::vector::const_iterator it = cluster.begin(), itend = cluster.end(); it < itend; ++it) { + // estimate start and end point + TCatPulseShape *pulse = *it; + + Double_t t = unit[i].Dot(pulse->GetPos() - centroid); + if (max < t) { + max = t; + itStartPoint = it; + } + if (t < min) { + min = t; + itEndPoint = it; + } + } - Double_t t = unit.Dot(pulse->GetPos() - centroid); + fProjectedLength[i] = TMath::Abs(unit[i].Dot((*itStartPoint)->GetPos() - (*itEndPoint)->GetPos())); - if (max < t) { - max = t; - itStartPoint = it; - } - if (t < min) { - min = t; - itEndPoint = it; + // take first principal axis is chosen for an estimation + if (i == 0) { + fStartPoint = (*itStartPoint)->GetPos(); + fEndPoint = (*itEndPoint)->GetPos(); } } - fStartPoint = centroid - ((*itStartPoint)->GetPos() - centroid).Mag() * unit; - fEndPoint = centroid + ((*itEndPoint)->GetPos() - centroid).Mag() * unit; - - fStartPoint = (*itStartPoint)->GetPos(); - fEndPoint = (*itEndPoint)->GetPos(); if (fStartPoint.X() * fPrincipalAxis.X() < 0) { // correction for end <-> start event TVector3 tmp = fEndPoint; @@ -217,6 +226,8 @@ void TCatTpcTrack::EstimateFromHits(const std::vector& clu Double_t vx = fPrincipalAxis.X(); Double_t vy = fPrincipalAxis.Y(); Double_t vz = fPrincipalAxis.Z(); + + for (std::vector::const_iterator it = cluster.begin(), itend = cluster.end(); it < itend; ++it) { TCatPulseShape *pulse = *it; const Double_t &c = pulse->GetCharge(); @@ -229,6 +240,11 @@ void TCatTpcTrack::EstimateFromHits(const std::vector& clu varx += c * dyz * dyz; vary += c * dzx * dzx; varz += c * dxy * dxy; + + // calculate length onto principal axis + + + } vary = TMath::Sqrt(vary/sum); varz = TMath::Sqrt(varz/sum); @@ -344,6 +360,8 @@ void TCatTpcTrack::Streamer(TBuffer &R__b) fDirection.Streamer(R__b); fPrincipalAxis.Streamer(R__b); fDeviation.Streamer(R__b); + if (R__v > 2) fProjectedLength.Streamer(R__b); + R__b >> fCharge; R__b >> fEnergyDeposit; R__b >> fRange; @@ -385,6 +403,7 @@ void TCatTpcTrack::Streamer(TBuffer &R__b) fDirection.Streamer(R__b); fPrincipalAxis.Streamer(R__b); fDeviation.Streamer(R__b); + fProjectedLength.Streamer(R__b); // for version 3 R__b << fCharge; R__b << fEnergyDeposit; R__b << fRange; diff --git a/TCatTpcTrack.h b/TCatTpcTrack.h index 2fb2daef8d61ebd3c763da15157483e667b25aa0..664dc6a63abd134ec190571b71f2c5e48bbbe788 100644 --- a/TCatTpcTrack.h +++ b/TCatTpcTrack.h @@ -3,7 +3,7 @@ * @brief container for tpc track * * @date Created : 2016-07-24 11:37:39 JST - * Last Modified : 2020-04-03 16:53:59 JST (ota) + * Last Modified : 2020-05-19 19:07:14 JST (ota) * @author Shinsuke OTA * * (C) 2016 Shinsuke OTA @@ -118,6 +118,8 @@ protected: TVector3 fPrincipalAxis; // principal axis of hits TVector3 fDeviation; // deviation of position from the principal axis + TVector3 fProjectedLength; // length projected onto each principal axis + Double_t fCharge; // sum of charge Double_t fEnergyDeposit; // sum of energy deposit @@ -139,7 +141,7 @@ protected: private: - ClassDef(TCatTpcTrack,2) // container for tracking result + ClassDef(TCatTpcTrack,3) // container for tracking result }; #endif // INCLUDE_GUARD_UUID_DF45729E_1BB8_49A2_A2ED_968FE7E7ED47