Functions | |
| Bool_t | IsLI (const NtpStRecord &ntp) |
| Bool_t | IsLIforNC (const NtpStRecord &ntp) |
|
|
This algorithm works by using known features of the topology of LI events. In particular it uses the asymmetry that comes from only flashing one side of the detector but also the crate boundaries. Definition at line 238 of file LISieve.cxx. References VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), RecHeader::GetVldContext(), MAXMSG, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRStrip::time0, and NtpSRStrip::time1. Referenced by EventQualAna::Analyze(), NuAnalysis::ChargeSignCut(), MadTVAnalysis::CreatePAN(), MadPIDAnalysis::CreatePAN(), NuAnalysis::EnergySpect(), NuExtraction::ExtractLITags(), ANtpInfoObjectFiller::FillHeaderInformation(), and NuAnalysis::LIRejectionTest(). 00239 {
00245
00246 //only analyse if FD
00247 if (ntp.GetHeader().GetVldContext().GetDetector()!=
00248 Detector::kFar) return false;
00249
00250 //variable to turn on all the useful messages if required
00251 Msg::LogLevel_t logLevel=Msg::kVerbose;
00252
00253 //this variable can turn off the noise removal in the event
00254 const Bool_t cleanTheEvent=true;
00255
00256 //store the set of times for the snarl
00257 multiset<Double_t> times;
00258
00259 Int_t numStrips=ntp.stp->GetEntries();
00260 if (numStrips<=0) return false;//pass the event, not LI
00262 //loop over the strips in the snarl
00264 for (Int_t hit=0;hit<numStrips;hit++){
00265 NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);
00266
00267 Double_t time=strip->time0;
00268 if (time<0 || time>1) time=strip->time1;
00269
00270 if (time>0 && time<=1) {
00271 times.insert(time);
00272 }
00273 //else just don't put the time in the map - clearly rubbish
00274 }//end of loop over strips
00275
00276 //get the median time from the map
00277 multiset<Double_t>::const_iterator it=times.begin();
00278 advance(it,times.size()/2);
00279 Double_t medianTime=*it;
00280 MAXMSG("LISieve",Msg::kDebug,100)
00281 <<"LISieve: Median time="<<medianTime<<endl;
00282
00283 map<Int_t,Int_t> hpp; // <plane, # of hits>
00284 Float_t sumSigCorEast=0;
00285 Float_t sumSigCorWest=0;
00286 //second loop
00287 for (Int_t hit=0;hit<numStrips;hit++){
00288 NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);
00289
00290 if (cleanTheEvent) {
00291 if (strip->ph0.sigcor > 0){
00292 if (strip->time0<medianTime-200e-9 ||
00293 strip->time0>medianTime+200e-9) continue;
00294 }
00295 if (strip->ph1.sigcor > 0){
00296 if (strip->time1<medianTime-200e-9 ||
00297 strip->time1>medianTime+200e-9) continue;
00298 }
00299 } // if (cleanTheEvent)
00300
00301 sumSigCorEast+=strip->ph0.sigcor;
00302 sumSigCorWest+=strip->ph1.sigcor;
00303
00304 hpp[strip->plane] += 1;
00305 }//end of loop over strips
00306
00307 Float_t avHpp=0;
00308 vector<Float_t> pb(8);//to store fractions hit in each pb region
00309
00310 //loop over all the planes flashed
00311 for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
00312 hppIt!=hpp.end();++hppIt){
00313 avHpp+=hppIt->second;
00314
00315 Int_t pl=hppIt->first;
00316
00317 if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
00318 else if (pl>=65 && pl<=128) pb[1]+=1./64;
00319 else if (pl>=129 && pl<=192) pb[2]+=1./64;
00320 else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
00321 else if (pl>=250 && pl<=313) pb[4]+=1./64;
00322 else if (pl>=314 && pl<=377) pb[5]+=1./64;
00323 else if (pl>=378 && pl<=441) pb[6]+=1./64;
00324 else if (pl>=442 && pl<=485) pb[7]+=1./44;
00325
00326 MAXMSG("LISieve",Msg::kVerbose,2000)
00327 <<"plane="<<pl<<", hpp="<<hppIt->second<<endl;
00328 }
00329 if (hpp.size()>0) avHpp/=hpp.size();
00330 Float_t sumSig=sumSigCorEast+sumSigCorWest;
00331 Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
00332 Float_t asym=0;
00333
00334 if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;
00335
00336 // Look for pulser box fractions based on planes only
00337
00338 Float_t fractFlashed1=0;
00339 Float_t fractFlashed2=0;
00340
00341 for (vector<Float_t>::iterator it=pb.begin();
00342 it!=pb.end();++it)
00343 {
00344 Float_t fract=*it;
00345 if (fract>fractFlashed1){
00346 //copy the 1 into 2
00347 fractFlashed2=fractFlashed1;
00348 //then set new highest fractFlashed
00349 fractFlashed1=fract;
00350 }
00351 //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
00352 else if (fract>fractFlashed2) {
00353 fractFlashed2=fract;
00354 }
00355
00356 MAXMSG("LISieve",Msg::kVerbose,200)
00357 <<"fract="<<fract<<", f1="<<fractFlashed1
00358 <<", f2="<<fractFlashed2<<endl;
00359 }
00360
00361 Float_t ratioFlashed=0;
00362 if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;
00363
00364
00365 Bool_t li=false;
00366
00367 // check if LI
00368 // don't check the asymmetry if the VA is saturating
00369 if ( (asym>0.5 || minSigCorEW>1.7e6) &&
00370 avHpp>3 && ratioFlashed<0.05 &&
00371 fractFlashed1>0.85 && fractFlashed2<0.1 )
00372 {
00373 li=true;
00374
00375 MAXMSG("LISieve",logLevel,100)
00376 <<"LI Event:"<<endl
00377 <<" Asymmetry="<<asym<<endl
00378 <<" Av. HPP ="<<avHpp<<endl
00379 <<" Fract Fl1="<<fractFlashed1
00380 <<", Fract Fl2="<<fractFlashed2
00381 <<" Ratio Fla="<<ratioFlashed << endl;
00382 }
00383 else
00384 {
00385 MAXMSG("LISieve",logLevel,10)
00386 <<"Non-LI Event:"<<endl
00387 <<" Asymmetry="<<asym<<endl
00388 <<" Av. HPP ="<<avHpp<<endl
00389 <<" Fract Fl1="<<fractFlashed1
00390 <<", Fract Fl2="<<fractFlashed2
00391 <<" Ratio Fla="<<ratioFlashed <<endl;
00392 }
00393
00394 return li;
00395 }
|
|
|
This algorithm works by using known features of the topology of LI events. In particular it uses the asymmetry that comes from only flashing one side of the detector but also the crate boundaries. Definition at line 22 of file LISieve.cxx. References VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), RecHeader::GetVldContext(), MAXMSG, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRPulseHeight::sigcor, NtpStRecord::stp, NtpSRStrip::time0, and NtpSRStrip::time1. Referenced by ANtpInfoObjectFiller::FillHeaderInformation(). 00023 {
00029
00030 //only analyse if FD
00031 if (ntp.GetHeader().GetVldContext().GetDetector()!=
00032 Detector::kFar) return false;
00033
00034 //variable to turn on all the useful messages if required
00035 Msg::LogLevel_t logLevel=Msg::kVerbose;
00036
00037 //this variable can turn off the noise removal in the event
00038 const Bool_t cleanTheEvent=true;
00039
00040 //store the set of times for the snarl
00041 multiset<Double_t> times;
00042
00043 Int_t numStrips=ntp.stp->GetEntries();
00044 if (numStrips<=0) return false;//pass the event, not LI
00046 //loop over the strips in the snarl
00048 for (Int_t hit=0;hit<numStrips;hit++){
00049 NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);
00050
00051 Double_t time=strip->time0;
00052 if (time<0 || time>1) time=strip->time1;
00053
00054 if (time>0 && time<=1) {
00055 times.insert(time);
00056 }
00057 //else just don't put the time in the map - clearly rubbish
00058 }//end of loop over strips
00059
00060 //get the median time from the map
00061 multiset<Double_t>::const_iterator it=times.begin();
00062 advance(it,times.size()/2);
00063 Double_t medianTime=*it;
00064 MAXMSG("LISieve",Msg::kDebug,100)
00065 <<"LISieve: Median time="<<medianTime<<endl;
00066
00067 map<Int_t,Int_t> hpp; // <plane, # of hits>
00068 Float_t sumSigCorEast=0;
00069 Float_t sumSigCorWest=0;
00070 //second loop
00071 for (Int_t hit=0;hit<numStrips;hit++){
00072 NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);
00073
00074 if (cleanTheEvent) {
00075 if (strip->ph0.sigcor > 0){
00076 if (strip->time0<medianTime-200e-9 ||
00077 strip->time0>medianTime+200e-9) continue;
00078 }
00079 if (strip->ph1.sigcor > 0){
00080 if (strip->time1<medianTime-200e-9 ||
00081 strip->time1>medianTime+200e-9) continue;
00082 }
00083 } // if (cleanTheEvent)
00084
00085 sumSigCorEast+=strip->ph0.sigcor;
00086 sumSigCorWest+=strip->ph1.sigcor;
00087
00088 hpp[strip->plane] += 1;
00089 }//end of loop over strips
00090
00091 Float_t avHpp=0;
00092 vector<Float_t> pb(8);//to store fractions hit in each pb region
00093 vector<Float_t> pbStrips(8);//to store fractionof strip hit in each pb region
00094
00095 //loop over all the planes flashed
00096 for (map<Int_t,Int_t>::iterator hppIt=hpp.begin();
00097 hppIt!=hpp.end();++hppIt){
00098 avHpp+=hppIt->second;
00099
00100 Int_t pl=hppIt->first;
00101
00102 if (pl>=1 && pl<=64) pb[0]+=1./64;//0=bookend
00103 else if (pl>=65 && pl<=128) pb[1]+=1./64;
00104 else if (pl>=129 && pl<=192) pb[2]+=1./64;
00105 else if (pl>=193 && pl<=248) pb[3]+=1./56;//249=bookend
00106 else if (pl>=250 && pl<=313) pb[4]+=1./64;
00107 else if (pl>=314 && pl<=377) pb[5]+=1./64;
00108 else if (pl>=378 && pl<=441) pb[6]+=1./64;
00109 else if (pl>=442 && pl<=485) pb[7]+=1./44;
00110
00111
00112 if (pl>=1 && pl<=64) pbStrips[0]+=avHpp*1./(64*192.);//0=bookend
00113 else if (pl>=65 && pl<=128) pbStrips[1]+=avHpp*1./(64*192.);
00114 else if (pl>=129 && pl<=192) pbStrips[2]+=avHpp*1./(64*192.);
00115 else if (pl>=193 && pl<=248) pbStrips[3]+=avHpp*1./(56*192.);//249=bookend
00116 else if (pl>=250 && pl<=313) pbStrips[4]+=avHpp*1./(64*192.);
00117 else if (pl>=314 && pl<=377) pbStrips[5]+=avHpp*1./(64*192.);
00118 else if (pl>=378 && pl<=441) pbStrips[6]+=avHpp*1./(64*192.);
00119 else if (pl>=442 && pl<=485) pbStrips[7]+=avHpp*1./(44*192.);
00120
00121
00122
00123 MAXMSG("LISieve",Msg::kVerbose,2000)
00124 <<"plane="<<pl<<", hpp="<<hppIt->second<<endl;
00125 }
00126 if (hpp.size()>0) avHpp/=hpp.size();
00127 Float_t sumSig=sumSigCorEast+sumSigCorWest;
00128 Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
00129 Float_t asym=0;
00130
00131 if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;
00132
00133 // Look for pulser box fractions based on planes only
00134
00135 Float_t fractFlashed1=0;
00136 Float_t fractFlashed2=0;
00137
00138 for (vector<Float_t>::iterator it=pb.begin();
00139 it!=pb.end();++it)
00140 {
00141 Float_t fract=*it;
00142 if (fract>fractFlashed1){
00143 //copy the 1 into 2
00144 fractFlashed2=fractFlashed1;
00145 //then set new highest fractFlashed
00146 fractFlashed1=fract;
00147 }
00148 //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
00149 else if (fract>fractFlashed2) {
00150 fractFlashed2=fract;
00151 }
00152
00153 MAXMSG("LISieve",Msg::kVerbose,200)
00154 <<"fract="<<fract<<", f1="<<fractFlashed1
00155 <<", f2="<<fractFlashed2<<endl;
00156 }
00157
00158 Float_t ratioFlashed=0;
00159 if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;
00160 // --------------------------
00161 // Same as above but for strip use fraction
00162 // added Nov 10,2009 by Jeff de Jong
00163 // Using strips offers abetter way to discriminate
00164 // LI from signal as LI lights up everything
00165
00166 Float_t fractStripsFlashed1=0;
00167 Float_t fractStripsFlashed2=0;
00168
00169 for (vector<Float_t>::iterator it=pbStrips.begin();
00170 it!=pbStrips.end();++it)
00171 {
00172 Float_t fract=*it;
00173 if (fract>fractStripsFlashed1){
00174 //copy the 1 into 2
00175 fractStripsFlashed2=fractStripsFlashed1;
00176 //then set new highest fractFlashed
00177 fractStripsFlashed1=fract;
00178 }
00179 //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
00180 else if (fract>fractStripsFlashed2) {
00181 fractStripsFlashed2=fract;
00182 }
00183
00184 MAXMSG("LISieve",Msg::kVerbose,200)
00185 <<"fractStrip="<<fract<<", f1="<<fractStripsFlashed1
00186 <<", f2="<<fractStripsFlashed2<<endl;
00187 }
00188
00189 Float_t ratioStripsFlashed=0;
00190 if (fractStripsFlashed1>0)
00191 ratioStripsFlashed=fractStripsFlashed2/fractStripsFlashed1;
00192
00193 //-------------------------
00194
00195
00196 Bool_t li=false;
00197
00198 // check if LI
00199 // don't check the asymmetry if the VA is saturating
00200 // if ( (asym>0.5 || minSigCorEW>1.7e6) &&
00201 // avHpp>3 && ratioFlashed<0.05 &&
00202 // fractFlashed1>0.85 && fractFlashed2<0.1 )
00203 // New cuts added on 10 Nov,2009 by J. de Jong
00204 if ( (asym>0.55 || minSigCorEW>1.7e6) &&
00205 ratioStripsFlashed<0.006 &&
00206 fractStripsFlashed1>0.02 )
00207 {
00208 li=true;
00209
00210 MAXMSG("LISieve",logLevel,100)
00211 <<"LI Event:"<<endl
00212 <<" Asymmetry="<<asym<<endl
00213 <<" Av. HPP ="<<avHpp<<endl
00214 <<" Fract Fl1="<<fractFlashed1
00215 <<", Fract Fl2="<<fractFlashed2
00216 <<" Ratio Fla="<<ratioFlashed
00217 <<" Strips Fl1="<<fractStripsFlashed1
00218 <<", StripsFl2="<<fractStripsFlashed2
00219 <<" Ratio Strips="<<ratioStripsFlashed<<endl;
00220 }
00221 else
00222 {
00223 MAXMSG("LISieve",logLevel,10)
00224 <<"Non-LI Event:"<<endl
00225 <<" Asymmetry="<<asym<<endl
00226 <<" Av. HPP ="<<avHpp<<endl
00227 <<" Fract Fl1="<<fractFlashed1
00228 <<", Fract Fl2="<<fractFlashed2
00229 <<" Ratio Fla="<<ratioFlashed
00230 <<" Strips Fl1="<<fractStripsFlashed1
00231 <<", StripsFl2="<<fractStripsFlashed2
00232 <<" Ratio Strips="<<ratioStripsFlashed<<endl;
00233 }
00234
00235 return li;
00236 }
|
1.3.9.1