loon: LISieve Namespace Reference

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.
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 
00071   //second loop
00072   for (Int_t hit=0;hit<numStrips;hit++){
00073     NtpSRStrip* strip = dynamic_cast<NtpSRStrip*>((*ntp.stp)[hit]);
00074 
00075     if (cleanTheEvent) {
00076       if (strip->ph0.sigcor > 0){
00077         if (strip->time0<medianTime-200e-9 || 
00078             strip->time0>medianTime+200e-9) continue;
00079       }
00080       if (strip->ph1.sigcor > 0){
00081         if (strip->time1<medianTime-200e-9 || 
00082             strip->time1>medianTime+200e-9) continue;
00083       }
00084     } // if (cleanTheEvent)
00085     
00086     sumSigCorEast+=strip->ph0.sigcor;
00087     sumSigCorWest+=strip->ph1.sigcor;
00088 
00089     hpp[strip->plane] += 1;
00090   }//end of loop over strips
00091 
00092   Float_t avHpp=0;
00093   vector<Float_t> pb(8);//to store fractions 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     MAXMSG("LISieve",Msg::kVerbose,2000)
00112       <<"plane="<<pl<<", hpp="<<hppIt->second<<endl;
00113   }
00114   if (hpp.size()>0) avHpp/=hpp.size();
00115   Float_t sumSig=sumSigCorEast+sumSigCorWest;
00116   Float_t minSigCorEW = TMath::Min(sumSigCorEast, sumSigCorWest);
00117   Float_t asym=0;
00118   
00119   if (sumSig) asym = (sumSig - 2 * minSigCorEW) / sumSig;
00120 
00121   Float_t fractFlashed1=0;
00122   Float_t fractFlashed2=0;
00123 
00124   for (vector<Float_t>::iterator it=pb.begin();
00125        it!=pb.end();++it){
00126     Float_t fract=*it;
00127     if (fract>fractFlashed1){
00128       //copy the 1 into 2
00129       fractFlashed2=fractFlashed1;
00130       //then set new highest fractFlashed
00131       fractFlashed1=fract;
00132     }
00133     //make sure to get fractFlashed2 if fractFlashed2<fract<frFla1
00134     else if (fract>fractFlashed2) {
00135       fractFlashed2=fract;
00136     }
00137 
00138     MAXMSG("LISieve",Msg::kVerbose,200)
00139       <<"fract="<<fract<<", f1="<<fractFlashed1
00140       <<", f2="<<fractFlashed2<<endl;
00141   }
00142 
00143   Float_t ratioFlashed=0;
00144   if (fractFlashed1>0) ratioFlashed=fractFlashed2/fractFlashed1;
00145 
00146   Bool_t li=false;
00147 
00148   // check if LI
00149   // don't check the asymmetry if the VA is saturating
00150   if ( (asym>0.5 || minSigCorEW>1.7e6) && 
00151        avHpp>3 && ratioFlashed<0.05 &&
00152        fractFlashed1>0.85 && fractFlashed2<0.1 ){
00153     li=true;
00154     
00155     MAXMSG("LISieve",logLevel,100)
00156       <<"LI Event:"<<endl
00157       <<"  Asymmetry="<<asym<<endl
00158       <<"  Av. HPP  ="<<avHpp<<endl
00159       <<"  Fract Fl1="<<fractFlashed1
00160       <<", Fract Fl2="<<fractFlashed2
00161       <<"  Ratio Fla="<<ratioFlashed<<endl;
00162   }
00163   else{
00164     MAXMSG("LISieve",logLevel,10)
00165       <<"Non-LI Event:"<<endl
00166       <<"  Asymmetry="<<asym<<endl
00167       <<"  Av. HPP  ="<<avHpp<<endl
00168       <<"  Fract Fl1="<<fractFlashed1
00169       <<", Fract Fl2="<<fractFlashed2
00170       <<"  Ratio Fla="<<ratioFlashed<<endl;
00171   }
00172   
00173   return li;
00174 }
LISieve Namespace Reference

Functions

Function Documentation


Functions
Bool_t	IsLI (const NtpStRecord &ntp)