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