#include <AlgChopListMitre.h>
Inheritance diagram for AlgChopListMitre:
Public Member Functions | |
| AlgChopListMitre () | |
| virtual | ~AlgChopListMitre () |
| virtual void | RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx) |
| virtual void | Trace (const char *c) const |
|
|
Definition at line 50 of file AlgChopListMitre.cxx. 00051 {
00052 }
|
|
|
Definition at line 55 of file AlgChopListMitre.cxx. 00056 {
00057 }
|
|
||||||||||||||||
|
Algorithm to chop by using the summed energy waveform of the whole calorimeter. Implements AlgBase. Definition at line 62 of file AlgChopListMitre.cxx. References CandHandle::AddDaughterLink(), digits(), RawRecord::FindRawBlock(), Form(), AlgFactory::GetAlgHandle(), CandContext::GetCandRecord(), CandContext::GetDataIn(), VldContext::GetDetector(), MomNavigator::GetFragment(), AlgFactory::GetInstance(), CandContext::GetMom(), Calibrator::GetTDCFromTime(), RecMinos::GetVldContext(), Calibrator::Instance(), k1pe, kQieRcid, PlexPlaneId::LastPlaneNearCalor(), CandDigitList::MakeCandidate(), MSG, and CandHandle::SetName(). 00065 {
00069
00070 assert(candHandle.InheritsFrom("CandChopListHandle"));
00071 CandChopListHandle &chopList = dynamic_cast<CandChopListHandle &>(candHandle);
00072
00073 assert(candContext.GetDataIn());
00074 // Check for CandDigitListHandle input
00075 if (!(candContext.GetDataIn()->InheritsFrom("CandDigitListHandle"))) {
00076 MSG("Chop",Msg::kWarning ) << "Data into AlgChopListSharp2 is not a digit list." << std::endl;
00077 }
00078
00079 const CandDigitListHandle *cdlh_ptr =
00080 dynamic_cast<const CandDigitListHandle*>(candContext.GetDataIn());
00081
00082 const MomNavigator *mom = candContext.GetMom();
00083 RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
00084 if (!rr) {
00085 MSG("Chop", Msg::kWarning) << "No RawRecord in MOM." << endl;
00086 return;
00087 }
00088 const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *>
00089 (rr->FindRawBlock("RawDigitDataBlock"));
00090 if (!rddb) {
00091 MSG("Chop", Msg::kWarning) << "No RawDigitDataBlock in RawRecord." << endl;
00092 return;
00093 }
00094
00095 // Get setup for the DigitList maker algorithm
00096 AlgFactory &af = AlgFactory::GetInstance();
00097 AlgHandle ah = af.GetAlgHandle("AlgChop","default");
00098 CandContext cxx(this,candContext.GetMom());
00099
00100 const VldContext &context = *(candContext.GetCandRecord()->GetVldContext());
00101 if(context.GetDetector() != DetectorType::kNear)
00102 MSG("Chop",Msg::kError) << "Running the Sharp2 algorithm on FD data is a no-no!" << endl;
00103
00104 Calibrator& cal = Calibrator::Instance();
00105 UgliGeomHandle ugli(context);
00106
00107 // Now do the actual slicing.
00108
00109 // First, make a nice stl vector of the digits.
00110 DigitVector digits(cdlh_ptr);
00111
00112 UInt_t ndigits = digits.size();
00113 UInt_t nchop = 0;
00114
00115 // Sort the list by time.
00116 // Not neccessary for this algorithm.
00117 // std::sort(digits.begin(), digits.end(), compareDigitTimes());
00118
00119
00120 // Also, I want some other pieces of info:
00121 std::vector<int> digit_tdc(ndigits);
00122 std::vector<UInt_t> digit_plane(ndigits);
00123 //std::vector<float> digit_tpos(ndigits);
00124 for(UInt_t i=0;i<ndigits;i++) {
00125 digit_tdc[i] = (cal.GetTDCFromTime(digits[i].GetTime(CalTimeType::kNone), kQieRcid));
00126 digit_plane[i] = digits[i].GetPlexSEIdAltL().GetPlane();
00127 //if(digit_plane[i]<=PlexPlaneId::LastPlaneNearCalor())
00128 // digit_tpos[i] = ugli.GetStripHandle(digits[i].GetPlexSEIdAltL().GetBestSEId()).GetTPos();
00129 //else
00130 // digit_tpos[i] = -999;
00131 }
00132
00133 // Find first and last times. Add some padding so sertain operations are easier to code.
00134 Int_t tfirst = digit_tdc[0];
00135 Int_t tlast = digit_tdc[0];
00136 for(UInt_t i=0;i<ndigits;i++) {
00137 if(digit_tdc[i] < tfirst) tfirst = digit_tdc[i];
00138 if(digit_tdc[i] > tlast ) tlast = digit_tdc[i];
00139 }
00140 tfirst-=5;
00141 tlast +=5;
00142
00143
00144 // Make the energy histogram.
00145 MSG("Chop",Msg::kDebug) << "Running Chop_Sharp2" << endl;
00146
00147 UInt_t numBins = tlast-tfirst;
00148
00149 // Create the energy-time profile(s).
00150
00151 std::vector< float > profAngles;
00152 profAngles.push_back(0);
00153 profAngles.push_back(1.0/60.); // 1 bucket per 60 planes.
00154 profAngles.push_back(-1.0/60.); // 1 bucket per 60 planes.
00155 UInt_t nProf = profAngles.size();
00156
00157 std::vector< std::vector<float> >profiles(nProf);
00158 std::vector<float> meanPlane(numBins,0);
00159
00160
00161 for(UInt_t p = 0; p<nProf; p++) {
00162 profiles[p].resize(numBins,0.);
00163 }
00164
00165 for(UInt_t idig = 0; idig < ndigits; idig++ ) {
00166 float sigcor = digits[idig].GetCharge(CalDigitType::kSigCorr);
00167 float tdcbin = digit_tdc[idig]-tfirst;
00168 float dplane = digit_plane[idig];
00169 if(dplane>PlexPlaneId::LastPlaneNearCalor())
00170 dplane = PlexPlaneId::LastPlaneNearCalor();
00171 dplane = dplane-60.;
00172
00173 std::vector<int> proj(nProf);
00174 for(UInt_t p = 0; p<nProf; p++)
00175 proj[p] = TMath::Nint(tdcbin + profAngles[p]*dplane);
00176
00177 if((tdcbin<0) || ((int)numBins<=tdcbin)) MSG("Chop",Msg::kDebug) << "Whups!" << endl;
00178 else if(digit_plane[idig]<=PlexPlaneId::LastPlaneNearCalor()) {
00179 meanPlane[proj[0]] += dplane*sigcor;
00180
00181 for(UInt_t p = 0; p<nProf; p++) {
00182 profiles[p][proj[p]] += sigcor;
00183 }
00184 }
00185 }
00186
00187 for(UInt_t i=0;i<numBins;i++) meanPlane[i] /= profiles[0][i];
00188
00189
00190 std::vector<int> cuts;
00191 std::vector<int> cuts_type;
00192
00193 // Insert a vertical cut at beginning of spill:
00194 cuts.push_back(0);
00195 cuts_type.push_back(0);
00196
00197
00198 int peak_bin = 0;
00199 bool rising = true;
00200 for(UInt_t bin=1;bin<numBins-1;bin++) {
00201 float dq = profiles[0][bin+1] - profiles[0][bin];
00202 float max_climb = 2.5*sqrt(fabs(profiles[0][bin])/k1pe)*k1pe;
00203 if(max_climb < 4.*k1pe) max_climb = max_climb*2.;
00204
00205 bool falling = !rising;
00206
00207 if(falling && ((bin-peak_bin)<8))
00208 if(max_climb<5.*k1pe) max_climb = 5.*k1pe;
00209
00210 if(falling) {
00211 if(dq > max_climb){
00212 cuts.push_back(bin);
00213 cuts_type.push_back(0);
00214 rising = true;
00215 }
00216 }
00217
00218 if(rising) {
00219 if(dq< 0){
00220 rising = false;
00221 peak_bin = bin;
00222 }
00223 }
00224 }
00225
00226 // We have first set of trial cuts. Now re-cut by the mitre saw.
00227 // Loop through cuts:
00228 for(UInt_t icut=0;icut<cuts.size(); icut++) {
00229 // find the lowest-charge bin of all the mitres, going back and forth +-1 bin
00230 float lowest_size = profiles[0][cuts[icut]];
00231 if(lowest_size>0) { // we can improve
00232
00233 MSG("Chop",Msg::kDebug) << "Evaluating cut " << icut
00234 << " at bin " << cuts[icut] << " mean plane " << meanPlane[cuts[icut]] << endl;
00235
00236 int a = cuts[icut]-5; if(a<0) a=0;
00237 int b = cuts[icut]+5; if(b>=(int)numBins) b=numBins-1;
00238 for(int i=a; i<=b; i++) {
00239 MSG("Chop",Msg::kDebug) << "Bin: " << i;
00240 for(UInt_t p = 0; p<nProf; p++) MSG("Chop",Msg::kDebug) << "\t" << Form("%6.0f",profiles[p][i]);
00241 MSG("Chop",Msg::kDebug) << endl;
00242 }
00243
00244 for(UInt_t p = 0; p<nProf; p++) {
00245 for(Int_t bin=cuts[icut]; bin<= cuts[icut]; bin++) {
00246 if(profiles[p][bin] < lowest_size) {
00247 lowest_size = profiles[p][bin];
00248 cuts[icut] = bin;
00249 cuts_type[icut] = p;
00250 }
00251 }
00252 }
00253
00254 MSG("Chop",Msg::kDebug) << " -> Cutting bin " << cuts[icut] << " on projection " << cuts_type[icut]
00255 << " lowbin: " << lowest_size << endl;
00256
00257 }
00258 }
00259
00260
00261 // Insert vertical cut at end:
00262 cuts.push_back(numBins-1);
00263 cuts_type.push_back(0);
00264
00265 for(UInt_t icut=0; icut<cuts.size()-1; icut++) {
00266 DigitVector chop;
00267
00268 for(UInt_t idig = 0; idig < ndigits; idig++ ) {
00269
00270 float tdcbin = digit_tdc[idig]-tfirst;
00271 float dplane = digit_plane[idig];
00272 if(dplane>PlexPlaneId::LastPlaneNearCalor())
00273 dplane = PlexPlaneId::LastPlaneNearCalor();
00274 dplane -= 60.;
00275
00276 std::vector<int> proj(nProf);
00277 for(UInt_t p = 0; p<nProf; p++)
00278 proj[p] = TMath::Nint(tdcbin + profAngles[p]*dplane);
00279
00280 // Is it after start cut?
00281 if(proj[cuts_type[icut]] >= cuts[icut]) {
00282
00283 // Is it after end cut?
00284 if(proj[cuts_type[icut+1]] < cuts[icut+1] ) {
00285
00286 // Add it to the chop.
00287 chop.push_back(digits[idig]);
00288
00289 }
00290 }
00291 }
00292
00293 cxx.SetDataIn(&(chop));
00294 CandDigitListHandle chopHandle = CandDigitList::MakeCandidate(ah,cxx);
00295 chopHandle.SetName(Form("Chop %d",nchop++));
00296 chopList.AddDaughterLink(chopHandle);
00297
00298 MSG("Chop",Msg::kDebug) << "Creating chop. "
00299 << " Start: " << cuts[icut] << " Stop: " << cuts[icut+1]
00300 << " Digits: " << chop.size()
00301 << endl;
00302
00303 }
00304 }
|
|
|
Reimplemented from AlgBase. Definition at line 307 of file AlgChopListMitre.cxx. 00308 {
00309 }
|
1.3.9.1