AlgCalDetPID Class Reference

#include <AlgCalDetPID.h>

Inheritance diagram for AlgCalDetPID:
AlgBase

List of all members.

Public Member Functions

 AlgCalDetPID ()
virtual ~AlgCalDetPID ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void Trace (const char *c) const

Private Member Functions

bool HitTimeAnal (const CandDigitListHandle *cdlh, Double_t toftime, Float_t lowwin, Float_t highwin, Int_t &nbefore, Int_t &nafter, Int_t &nout, Int_t &nin, Float_t &chi2) const
bool IsAnOverlap (AlgConfig &ac, Int_t nbefore, Int_t nafter, Int_t nout, Int_t nin)
bool CERInTimeWin (AlgConfig &ac, const CandCalDetSIHandle *csih, Double_t abstime, Double_t toftime, Float_t win0low, Float_t win0high, Float_t win1low, Float_t win1high, Float_t win2low, Float_t win2high, bool &cer0, bool &cer1, bool &cer2)
Float_t OverlapLikelihood (TH1F *ehfar, TH1F *ehnear) const

Private Attributes

Bool_t fFirstEvent
Int_t fNAllowedBefore
Int_t fNAllowedAfter
Int_t fNAllowedOut
Int_t fNRequiredIn
Int_t fTimeFileNum
TFile * fTimeFile
TH1F * fThist [2]

Detailed Description

Definition at line 19 of file AlgCalDetPID.h.


Constructor & Destructor Documentation

AlgCalDetPID::AlgCalDetPID (  ) 

Definition at line 46 of file AlgCalDetPID.cxx.

00046                           :
00047    fFirstEvent(kTRUE),
00048    fNAllowedBefore(1), 
00049    fNAllowedAfter(1),
00050    fNAllowedOut(2), 
00051    fNRequiredIn(1), 
00052    fTimeFileNum(0)
00053 {
00054    fTimeFile=0;
00055    fThist[0]=0;
00056    fThist[1]=0;
00057 
00058 }

AlgCalDetPID::~AlgCalDetPID (  )  [virtual]

Definition at line 61 of file AlgCalDetPID.cxx.

References fTimeFile.

00062 {
00063    if(fTimeFile){
00064       if(fTimeFile->IsOpen()){
00065          fTimeFile->Close();
00066          delete fTimeFile;
00067          fTimeFile = 0;
00068       }
00069    }
00070 
00071 }


Member Function Documentation

bool AlgCalDetPID::CERInTimeWin ( AlgConfig ac,
const CandCalDetSIHandle csih,
Double_t  abstime,
Double_t  toftime,
Float_t  win0low,
Float_t  win0high,
Float_t  win1low,
Float_t  win1high,
Float_t  win2low,
Float_t  win2high,
bool &  cer0,
bool &  cer1,
bool &  cer2 
) [private]

Definition at line 601 of file AlgCalDetPID.cxx.

References CandCalDetSIHandle::GetKovADC1(), CandCalDetSIHandle::GetKovADC2(), CandCalDetSIHandle::GetKovADC3(), CandCalDetSIHandle::GetKovTimeStamp1(), CandCalDetSIHandle::GetKovTimeStamp2(), CandCalDetSIHandle::GetKovTimeStamp3(), Msg::kDebug, Msg::kError, MSG, Munits::ns, and tdc_convert.

Referenced by RunAlg().

00607                                                                    {
00608      // returns in cerX true if counter X passes the selection
00609      // returns cer0&cer1&cer2
00610 
00611      cer0=cer1=cer2=false;
00612      if(csih==0){
00613           MSG("CalDetPID",Msg::kError)
00614                <<"CERInTimeWin got a null CandCalDetSIHandle!"
00615                <<"Returning false!"<<endl;
00616           return false;
00617      }
00618      else{
00619           // in principle could use algorithm config to avoid considering
00620           // one of the counters for some particles
00621 
00622           // candcaldetsihandle starts counting cer counters at 1
00623           const Double_t tdc_convert=1.5625;
00624 
00625 
00626           Int_t cer0adc=csih->GetKovADC2(); // yes, I really mean this
00627           Float_t cer0time= toftime-
00628               (csih->GetKovTimeStamp2()*tdc_convert -abstime/Munits::ns); // too
00629           Int_t cer1adc=csih->GetKovADC1(); // too
00630           Float_t cer1time= toftime - 
00631               (csih->GetKovTimeStamp1()*tdc_convert -abstime/Munits::ns);//too
00632           Int_t cer2adc=csih->GetKovADC3();
00633           Float_t cer2time= toftime - 
00634                (csih->GetKovTimeStamp3()*tdc_convert-abstime/Munits::ns);
00635 
00636           MSG("CalDetPID",Msg::kDebug)<<"toftime - (kov2*tdc_convert -abstime/Munits::ns)= "
00637                                          <<toftime<<"- ("
00638                                          <<csih->GetKovTimeStamp2()*tdc_convert
00639                                          <<" - "<<abstime/Munits::ns<<")"
00640                                          <<" = "<<toftime -(csih->GetKovTimeStamp2()*tdc_convert -abstime/Munits::ns)
00641                                          <<endl;
00642           
00643           MSG("CalDetPID",Msg::kDebug)<<"toftime - (kov1*tdc_convert -abstime/Munits::ns)= "
00644                                          <<toftime<<"- ("
00645                                          <<csih->GetKovTimeStamp1()*tdc_convert
00646                                          <<" - "<<abstime/Munits::ns<<")"
00647                                          <<" = "<<toftime -(csih->GetKovTimeStamp1()*tdc_convert -abstime/Munits::ns)
00648                                          <<endl;
00649 
00650           MSG("CalDetPID",Msg::kDebug)<<"toftime - (kov3*tdc_convert -abstime/Munits::ns)= "
00651                                          <<toftime<<"- ("
00652                                          <<csih->GetKovTimeStamp3()*tdc_convert
00653                                          <<" - "<<abstime/Munits::ns<<")"
00654                                          <<" = "<<toftime -(csih->GetKovTimeStamp3()*tdc_convert -abstime/Munits::ns)
00655                                          <<endl;
00656           const float nullwin=100000;
00657           // set cerX true if counter X passes the selection
00658           if(cer0adc>0){
00659               if(win0low==-nullwin && win0high==nullwin) cer0=true;
00660               else if((cer0time<win0low) || (cer0time>win0high)) cer0=false;
00661               else cer0=true;
00662           }
00663           else cer0=true;
00664 
00665           if(cer1adc>0){
00666               if(win1low==-nullwin && win1high==nullwin) cer1=true;
00667               else if((cer1time<win1low) || (cer1time>win1high)) cer1=false;
00668               else cer1=true;
00669           }
00670           else cer1=true;
00671 
00672           if(cer2adc>0){
00673               if(win2low==-nullwin && win2high==nullwin) cer2=true;
00674               else if((cer2time<win2low) || (cer2time>win2high)) cer2=false;
00675               else cer2=true;
00676           }
00677           else cer2=true;
00678           
00679           return (cer0&cer1&cer2);
00680      }    
00681 
00682 }

bool AlgCalDetPID::HitTimeAnal ( const CandDigitListHandle cdlh,
Double_t  toftime,
Float_t  lowwin,
Float_t  highwin,
Int_t &  nbefore,
Int_t &  nafter,
Int_t &  nout,
Int_t &  nin,
Float_t &  chi2 
) const [private]

Definition at line 487 of file AlgCalDetPID.cxx.

References fThist, CandHandle::GetDaughterIterator(), PlexHandle::GetReadoutType(), Msg::kDebug, Msg::kError, ElecType::kQIE, ReadoutType::kScintStrip, CalTimeType::kT0, PlaneView::kU, PlaneView::kV, ElecType::kVA, MSG, Munits::nanosecond, and OverlapLikelihood().

Referenced by RunAlg().

00492 {
00493           // returns false if there is an error
00494           
00495           // zero counters
00496           nbefore=0; nafter=0; nout=0; nin=0;
00497           
00498           //make histograms to hold timing of current event
00499           int NBFAR = fThist[0]->GetNbinsX();
00500           Axis_t LOWFAR = fThist[0]->GetBinLowEdge(1);
00501           Axis_t HIFAR = fThist[0]->GetBinLowEdge(NBFAR)+
00502              fThist[0]->GetBinWidth(NBFAR);
00503 
00504           int NBNEAR = fThist[1]->GetNbinsX();
00505           Axis_t LOWNEAR = fThist[1]->GetBinLowEdge(1);
00506           Axis_t HINEAR = fThist[1]->GetBinLowEdge(NBNEAR)+
00507              fThist[1]->GetBinWidth(NBNEAR);
00508 
00509           MSG("CalDetPID",Msg::kDebug)<<"making event thists "
00510                                       <<" NBFAR "<<NBFAR
00511                                       <<" LOWFAR "<<LOWFAR
00512                                       <<" HIFAR "<<HIFAR
00513                                       <<" NBNEAR "<<NBNEAR
00514                                       <<" LOWNEAR "<<LOWNEAR
00515                                       <<" HINEAR "<<HINEAR<<endl;
00516 
00517           TH1F ehfar("ehfar","ehfar",NBFAR,LOWFAR,HIFAR);
00518           TH1F ehnear("ehnear","ehnear",NBNEAR,LOWNEAR,HINEAR);
00519   
00520           // check digit list
00521           CandDigitHandleItr hi(cdlh->GetDaughterIterator());
00522           if(!hi.IsValid()){
00523                MSG("CalDetPID",Msg::kError)<<"DaugherIter not valid"<<endl;
00524                return false;    
00525           }
00526           //begin looping over hits
00527           for(hi.Reset();hi.IsValid();hi.Next()){
00528               // to check that it isn't a cosmic counter or unknown
00529               const PlaneView::PlaneView_t pv=
00530                    (*hi)->GetPlexSEIdAltL().GetPlaneView();
00531               // to check that it IS connected to scint
00532               const PlexHandle ph(*((*hi)->GetVldContext()));
00533               const ReadoutType::Readout_t rt= 
00534                    ph.GetReadoutType((*hi)->GetChannelId());
00535               const ElecType::Elec_t et = (*hi)->GetChannelId().GetElecType();
00536 /*
00537   // the logic used to be
00538               if((pv==PlaneView::kA)||(pv==PlaneView::kB)
00539                  ||(pv==PlaneView::kUnknown)) {
00540 */
00541 /*  changing the logic so that it will consider near detector hits
00542               if( !(pv==PlaneView::kU || pv==PlaneView::kV) 
00543                   || rt!=ReadoutType::kScintStrip || (et!=ElecType::kVA)){
00544                   // don't do anything with this hit
00545                    continue;
00546               }
00547 */
00548               if( !(pv==PlaneView::kU || pv==PlaneView::kV) 
00549                   || rt!=ReadoutType::kScintStrip){
00550                   // don't do anything with this hit
00551                    continue;
00552               }
00553 
00554               Float_t time = (*hi)->GetSubtractedTime(CalTimeType::kT0);
00555               time/=Munits::nanosecond;
00556               MSG("CalDetPID",Msg::kDebug)<<"hittime - toftime= "
00557                                              <<time<<"-"<<toftime
00558                                              <<" = "<<time-toftime<<endl;
00559               time-=toftime;
00560               if(et==ElecType::kVA){
00561                  ehfar.Fill(time);
00562                  if(time<lowwin){nbefore++; nout++;}
00563                  else if(time>highwin){nafter++; nout++;}
00564                  else nin++;
00565               }
00566               else if(et==ElecType::kQIE){
00567                  ehnear.Fill(time);
00568               }
00569           }
00570           chi2 = OverlapLikelihood(&ehfar,&ehnear);
00571           return true;
00572 }

bool AlgCalDetPID::IsAnOverlap ( AlgConfig ac,
Int_t  nbefore,
Int_t  nafter,
Int_t  nout,
Int_t  nin 
) [private]

Definition at line 574 of file AlgCalDetPID.cxx.

References fNAllowedAfter, fNAllowedBefore, fNAllowedOut, fNRequiredIn, and Registry::Get().

Referenced by RunAlg().

00577 {
00578      // returns false if there is no overlap
00579      // returns true if there is an overlap
00580      // get the number of hits allowed out of the time interval from
00581      // AlgConfig
00582      
00583      Int_t tmpi;
00584      if(ac.Get("NAllowedBefore", tmpi)) fNAllowedBefore=tmpi;
00585      if(ac.Get("NAllowedAfter", tmpi)) fNAllowedAfter=tmpi;
00586      if(ac.Get("NAllowedOut", tmpi)) fNAllowedOut=tmpi;
00587      if(ac.Get("NRequiredIn", tmpi)) fNRequiredIn=tmpi;
00588    
00589      if(nout<=fNAllowedOut&&
00590         nbefore<=fNAllowedBefore&&
00591         nafter<=fNAllowedAfter&&
00592         nin>=fNRequiredIn){
00593           return false;
00594      }
00595      else{
00596           return true;
00597      }  
00598 
00599 }

Float_t AlgCalDetPID::OverlapLikelihood ( TH1F *  ehfar,
TH1F *  ehnear 
) const [private]

Definition at line 684 of file AlgCalDetPID.cxx.

References fThist, Msg::kDebug, and MSG.

Referenced by HitTimeAnal().

00685 {
00686    //now loop over the histogram bins to calculate the chi2
00687    Int_t nhits=0;
00688    Float_t flsum=0.;
00689    
00690    TH1F *eventhist[2];
00691    eventhist[0] = ehfar;
00692    eventhist[1] = ehnear;
00693 
00694    for(int i=0;i<2;i++){
00695       const int nevents = (const int)eventhist[i]->Integral();
00696       nhits+=nevents;
00697       for(int j=1;j<=eventhist[i]->GetNbinsX();j++){
00698          if(eventhist[i]->GetBinContent(j)>0){
00699             int x = (int)eventhist[i]->GetBinContent(j);
00700             float mean = fThist[i]->GetBinContent(j)*nevents;
00701             MSG("CalDetPID",Msg::kDebug)<<" crate "<<i
00702                                    <<" at time "<<eventhist[i]->GetBinCenter(j)
00703                                    <<" x "<<x<<" mean "<<mean<<endl;
00704 
00705             if(mean==0){
00706                continue;
00707             }
00708 
00709             double lnprob =  1.*x*log(mean)-mean;
00710             double lnnorm = 1.*x*log(1.*x)-1.*x;
00711 
00712             flsum+=(lnprob-lnnorm);
00713             
00714             MSG("CalDetPID",Msg::kDebug)<<"i "<<i
00715                                    <<" ln prob "<<lnprob
00716                                    <<" ln norm "<<lnnorm
00717                                    <<" fl sum "<<flsum<<endl;
00718          }
00719       }
00720       MSG("CalDetPID",Msg::kDebug)<<"i "<<i
00721                              <<" flsum "<<flsum<<endl;
00722    }
00723    if(nhits>0){
00724       flsum = -2*(flsum/(1.*nhits));  
00725    }
00726    else{
00727       flsum = 0;
00728    }
00729    return flsum;
00730 }

void AlgCalDetPID::RunAlg ( AlgConfig ac,
CandHandle ch,
CandContext cx 
) [virtual]

Implements AlgBase.

Definition at line 74 of file AlgCalDetPID.cxx.

References bfld::AsString(), CERInTimeWin(), fFirstEvent, RecDataRecord< T >::FindComponent(), fname, fThist, fTimeFile, fTimeFileNum, Registry::Get(), CandDigitListHandle::GetAbsTime(), CandContext::GetCandRecord(), CalDetCERRange::GetCER0Low(), CalDetCERRange::GetCER1Low(), CalDetCERRange::GetCER2Low(), CandContext::GetDataIn(), MomNavigator::GetFragment(), CandCalDetSIHandle::GetKovADC1(), CandCalDetSIHandle::GetKovADC2(), CandCalDetSIHandle::GetKovADC3(), CandContext::GetMom(), DbiResultPtr< T >::GetNumRows(), CalDetOverlapWin::GetParticleType(), CalDetCERTimeWin::GetParticleType(), CalDetCERRange::GetParticleType(), CalDetTOFRange::GetParticleType(), DbiResultPtr< T >::GetRow(), VldContext::GetSimFlag(), CalDetTOFRange::GetTDC0Low(), CalDetTOFRange::GetTDC2Low(), CalDetTOFRange::GetTDC2MinusTDC0Low(), CandCalDetSIHandle::GetTofTDC0(), CandCalDetSIHandle::GetTofTDC1(), CandCalDetSIHandle::GetTofTDC2(), CandCalDetSIHandle::GetTofTimeStamp(), DbiResultPtr< T >::GetValidityRec(), CalDetCERTimeWin::GetWin0High(), CalDetCERTimeWin::GetWin0Low(), CalDetCERTimeWin::GetWin1High(), CalDetCERTimeWin::GetWin1Low(), CalDetCERTimeWin::GetWin2High(), CalDetCERTimeWin::GetWin2Low(), CalDetOverlapWin::GetWinHigh(), CalDetOverlapWin::GetWinLow(), HitTimeAnal(), IsAnOverlap(), SimFlag::kData, Msg::kDebug, CalDetParticleType::kElectron, Msg::kError, Msg::kFatal, Msg::kInfo, CalDetParticleType::kKaon, SimFlag::kMC, CalDetParticleType::kMuon, CalDetParticleType::kPion, CalDetParticleType::kProton, SimFlag::kReroot, CalDetParticleType::kUnknown, Msg::kWarning, MSG, and Munits::ns.

00076 {
00077      bool ok=true;
00078 
00079    MSG("CalDetPID", Msg::kDebug) 
00080         << "Starting AlgCalDetPID::RunAlg()" << endl;
00081 
00082    if(fFirstEvent){
00083       //get filename from alg config
00084       Int_t tfn;
00085       if(ac.Get("TimeFileNum", tfn)){
00086          fTimeFileNum = tfn;
00087          MSG("CalDetPID", Msg::kDebug)<<"got TimeFileNum from config "
00088                                       <<fTimeFileNum<<endl;
00089       }
00090       else{
00091          fTimeFileNum = 0;
00092          MSG("CalDetPID", Msg::kError)<<"Could not get TimeFileName from config."
00093                                       <<"  Using default "<<fTimeFileNum
00094                                       <<" instead"<<endl;
00095       }
00096       string fTimeFileName;
00097       switch(fTimeFileNum){
00098       case 0: //far-far readout t11 2002
00099          fTimeFileName = "timinghists-2002-t11.root";
00100          break;
00101       case 1: //far-far readout t7 2002
00102          fTimeFileName = "timinghists-2002-t7.root";
00103          break;
00104       case 2: //far-near readout 2003
00105          fTimeFileName = "timinghists-2003-t7-nf.root";
00106          break;
00107       case 3: //near readout t7 2003
00108          fTimeFileName = "timinghists-2003-t7-n.root";
00109          break;
00110       case 4: //near readout t11 2003
00111          fTimeFileName = "timinghists-2003-t11.root";
00112          break;
00113       default:
00114          fTimeFileName = "timinghists.root";
00115          break;
00116       }
00117       
00118       // open input root file
00119       string fname(getenv("SRT_PRIVATE_CONTEXT"));
00120       fname+=("/CalDetPID/data/"+fTimeFileName);
00121       
00122       fTimeFile = new TFile(fname.c_str());
00123       if(!fTimeFile->IsOpen()){
00124          MSG("CalDetPID",Msg::kInfo)<<"Could not find "<<fname<<" in the base release. Will try the test release."<<endl;
00125          string fname_save=fname;
00126          fname=getenv("SRT_PUBLIC_CONTEXT");
00127          fname+="/CalDetPID/data/"+fTimeFileName;
00128          //crashes if it get's here because now fTimeFile goes out of scope
00129          //tv 6-16-2004
00130 //       TFile *fTimeFile = new TFile(fname.c_str());
00131          fTimeFile = new TFile(fname.c_str());
00132          if(!fTimeFile->IsOpen()){
00133             MSG("CalDetPID",Msg::kError)<<"Could not open timing hist file"
00134                                         <<fname<<endl
00135                                         <<"Also tried "<<fname_save<<endl;
00136             return;
00137          }
00138       }
00139       MSG("CalDetPID",Msg::kInfo)<<"Opened timing hist file "<<fname<<endl;
00140       //      fTimeFile->Print();
00141       for(int i=0;i<2;i++){
00142          char hn[50];
00143          sprintf(hn,"all_thist_c%d",i);
00144          fThist[i] = (TH1F *)fTimeFile->Get(hn);
00145          if(fThist[i]!=0){
00146             MSG("CalDetPID",Msg::kDebug)<<"Found histogram "<<hn<<endl;
00147          }
00148          else{
00149             MSG("CalDetPID",Msg::kError)<<"Could not find histogram "<<hn<<endl;
00150          }
00151       }
00152       fFirstEvent=kFALSE;
00153    }
00154 
00155 
00156    assert(ch.InheritsFrom("CandCalDetPIDHandle"));
00157    CandCalDetPIDHandle& cpidh = static_cast<CandCalDetPIDHandle&>(ch);
00158 
00159    assert(cx.GetDataIn());
00160    const CandRecord* crec = cx.GetCandRecord();
00161    // Get VldContext
00162    const VldContext vc = *(crec->GetVldContext());
00163    // is this data?
00164    const bool isdata = (vc.GetSimFlag()&SimFlag::kData) ? true : false;
00165    // is this mc?
00166    const bool ismc = ((vc.GetSimFlag()&SimFlag::kMC) ||
00167                       (vc.GetSimFlag()&SimFlag::kReroot)) ? true : false;
00168    
00169    if(isdata){
00170         // for data we must have a caldetsi candidate
00171         const CandCalDetSIHandle *csih =
00172              dynamic_cast<const CandCalDetSIHandle *> (cx.GetDataIn());
00173         assert(csih);
00174         // we must also have some digits
00175         const CandDigitListHandle* cdlh=
00176              dynamic_cast<const CandDigitListHandle*>
00177              (crec->FindCandHandle("CandDigitListHandle","canddigitlist"));
00178         if(cdlh == 0) {//failure
00179              MSG("CalDetPID",Msg::kFatal)<<"cdlh==0"<<endl;
00180              assert(0);
00181         }
00182 
00183         const Double_t toftime=
00184              (csih->GetTofTimeStamp()*1.5625 - cdlh->GetAbsTime()/Munits::ns);
00185 
00186         // talk to the db
00187    
00188         DbiResultPtr<CalDetTOFRange> rp_tof(vc,0);
00189         if(!rp_tof.GetValidityRec()){
00190              MSG("CalDetPID", Msg::kError)
00191                   <<"Could not look up selection ranges "
00192                   <<"from CalDetTOFRange table!\n"
00193                   <<"Bad DbiResultPtr!?"
00194                   <<endl;
00195              ok=false;
00196         }
00197 
00198         DbiResultPtr<CalDetCERRange> rp_cer(vc,0);
00199         if(!rp_cer.GetValidityRec()){
00200              MSG("CalDetPID", Msg::kError)
00201              <<"Could not look up selection ranges "
00202              <<"from CalDetCERRange table!\n"
00203              <<"Bad DbiResultPtr!?"
00204              <<endl;
00205              ok=false;
00206         }
00207 
00208 
00209    
00210         // The algorithm:
00211         // Look through each entry in the db.  
00212         // There should not be many of them.
00213         // One entry may correspond to one or more types of particles.
00214         // The PARTICLETYPE entry in the table rows will be used to set high
00215         // the corresponding bits in case the data is in the appropriate range.
00216         // This can be done seperately for tof and cer with the results anded
00217         // at the end.
00218 
00219         // do id by tof
00220         UInt_t tof_pid=0;       
00221         if(ok){   
00222              // here is the tof data
00223              const Int_t tdc0=csih->GetTofTDC0();
00224              const Int_t tdc1=csih->GetTofTDC1();
00225              const Int_t tdc2=csih->GetTofTDC2();
00226              MSG("CalDetPID", Msg::kDebug)
00227                   <<"TOF 0 1 2 = "<<tdc0<<" "<<tdc1<<" "<<tdc2<<endl;
00228              // loop over rows
00229              const int nrows_tof = rp_tof.GetNumRows();
00230              for(int i=0; i<nrows_tof; i++){
00231                   // get a row
00232                   const CalDetTOFRange* tr = rp_tof.GetRow(i);
00233                   if(!tr){
00234                        // something is really hosed
00235                        MSG("CalDetPID", Msg::kFatal)
00236                             <<"Could not look up row "<<i
00237                             <<" in CalDetTOFRange table!"
00238                             <<endl;
00239                        ok=false;
00240                        assert(0);
00241                   }
00242                   else{
00243                        // for now, only test on tdc2 and tdc0
00244                        // tdc1 was never used...
00245                        // when we ran middle and ds counters the middle
00246                        // was plugged into tdc2
00247                        if((tdc0 < tr->GetTDC0High()) 
00248                           && (tdc0 > tr->GetTDC0Low())
00249                           && (tdc2 < tr->GetTDC2High()) 
00250                           && (tdc2 > tr->GetTDC2Low())
00251                           && (tdc2-tdc0 < tr->GetTDC2MinusTDC0High())
00252                           && (tdc2-tdc0 > tr->GetTDC2MinusTDC0Low())
00253                             )
00254                             tof_pid|=(tr->GetParticleType());
00255                   }
00256              }
00257         }
00258 
00259         // do id by cer
00260         UInt_t cer_pid=0;
00261         if(ok){
00262              // here is the cer data
00263              const Int_t cer0=csih->GetKovADC2();
00264              const Int_t cer1=csih->GetKovADC1();
00265              const Int_t cer2=csih->GetKovADC3();
00266              MSG("CalDetPID", Msg::kDebug)
00267                   <<"CER 0 1 2 = "<<cer0<<" "<<cer1<<" "<<cer2<<endl;
00268              // loop over rows
00269              const int nrows_cer = rp_cer.GetNumRows();
00270              for(int i=0; i<nrows_cer; i++){
00271                   const CalDetCERRange* cr = rp_cer.GetRow(i);
00272                   if(!cr){
00273                        // something is really hosed
00274                        MSG("CalDetPID", Msg::kFatal)
00275                             <<"Could not look up row "<<i
00276                             <<" in  CalDetCERRange table!"
00277                             <<endl;
00278                        ok=false;
00279                        assert(0);
00280                   }
00281                   else{  
00282                        if((cer0<cr->GetCER0High())
00283                           &&(cer0>cr->GetCER0Low())
00284                           &&(cer1<cr->GetCER1High())
00285                           &&(cer1>cr->GetCER1Low())
00286                           &&(cer2<cr->GetCER2High())
00287                           &&(cer2>cr->GetCER2Low())
00288                             ) 
00289                             cer_pid|=(cr->GetParticleType());
00290                   }
00291              }    
00292         }
00293 
00294         MSG("CalDetPID", Msg::kDebug)
00295              <<"tof_pid = "<<hex<<tof_pid<<dec
00296              <<",  Identified via TOF as "
00297              <<CalDetParticleType::AsString(tof_pid)<<endl;
00298         
00299         MSG("CalDetPID", Msg::kDebug)
00300              <<"cer_pid = "<<hex<<cer_pid<<dec
00301              <<",  Identified via CER as "
00302              <<CalDetParticleType::AsString(cer_pid)<<endl;
00303         
00304         // Now we have the basic pid information.
00305         // Look to see if this is an overlapper.
00306 
00307         // For things that are not overlappers, set the appropriate bits high
00308         // otherwise do not set them
00309         // If we cannot talk to the db then cannot_run_overlap=true;
00310         UInt_t overlap_pid=0;
00311         bool cannot_run_overlap=false;
00312         Float_t chi2 = -1.;
00313         DbiResultPtr<CalDetOverlapWin> rp_cow(vc,0);
00314         if(!rp_cow.GetValidityRec()){
00315              MSG("CalDetPID", Msg::kWarning)
00316                   <<"Could not look up selection ranges "
00317                   <<"from CalDetOverlapWin table!\n"
00318                   <<"Bad DbiResultPtr!?"
00319                   <<endl;
00320              // but maybe it could live with this?
00321              cannot_run_overlap=true;
00322              ok=false;
00323         }
00324         else{
00325              const int nrows_cow=rp_cow.GetNumRows();
00326              for(int i=0; i<nrows_cow; i++){
00327                   const CalDetOverlapWin* cow = rp_cow.GetRow(i);
00328                   if(!cow){
00329                        // something is really hosed
00330                        MSG("CalDetPID", Msg::kFatal)
00331                             <<"Could not look up row "<<i
00332                             <<" in  CalDetOverlapWin table!"
00333                             <<endl;
00334                        ok=false;
00335                        assert(0);
00336                   }
00337                   else{
00338                        const Float_t lowwin = cow->GetWinLow();
00339                        const Float_t highwin = cow->GetWinHigh();
00340                        Int_t nbefore,nafter,nout,nin;
00341                        HitTimeAnal(cdlh,toftime,lowwin,highwin,
00342                                    nbefore,nafter,nout,nin,chi2);
00343                        bool overlap = IsAnOverlap(ac,nbefore,nafter,
00344                                                   nout,nin);
00345                        if(!overlap){
00346                             overlap_pid|=(cow->GetParticleType());
00347                        }
00348                   }
00349              }
00350         }
00351                 
00352         //
00353         // now work on certimewin
00354         UInt_t ctw_pid=0;
00355         bool cannot_run_ctw=false;
00356         DbiResultPtr<CalDetCERTimeWin> rp_ctw(vc,0);
00357         if(!rp_ctw.GetValidityRec()){
00358              MSG("CalDetPID", Msg::kError)
00359                   <<"Could not look up selection ranges "
00360                   <<"from CalDetCERTimeWin table!\n"
00361                   <<"Bad DbiResultPtr!?"
00362                   <<endl;
00363              // but maybe it could live with this?
00364              cannot_run_ctw=true;
00365              ok=false;
00366         }
00367         else{
00368              const int nrows_ctw=rp_ctw.GetNumRows();
00369              for(int i=0; i<nrows_ctw; i++){
00370                   const CalDetCERTimeWin* ctw = rp_ctw.GetRow(i);
00371                   if(!ctw){
00372                        // something is really hosed
00373                        MSG("CalDetPID", Msg::kFatal)
00374                             <<"Could not look up row "<<i
00375                             <<" in  CalDetCERTimeWin table!"
00376                             <<endl;
00377                        ok=false;
00378                        assert(0);
00379                   }
00380                   else{
00381                        const Float_t win0low = ctw->GetWin0Low();
00382                        const Float_t win0high = ctw->GetWin0High();
00383                        const Float_t win1low = ctw->GetWin1Low();
00384                        const Float_t win1high = ctw->GetWin1High();
00385                        const Float_t win2low = ctw->GetWin2Low();
00386                        const Float_t win2high = ctw->GetWin2High();
00387                        bool cer0=false;
00388                        bool cer1=false;
00389                        bool cer2=false;
00390                        MSG("CalDetPID", Msg::kDebug)
00391                            <<(*ctw)<<endl;
00392                        const bool in_win = CERInTimeWin(ac, csih,
00393                                                         cdlh->GetAbsTime(),
00394                                                         toftime,
00395                                                         win0low,win0high,
00396                                                         win1low,win1high,
00397                                                         win2low,win2high,
00398                                                         cer0, cer1, cer2);
00399                        if(in_win){
00400                             ctw_pid|=(ctw->GetParticleType());
00401                        }
00402                   }
00403              }
00404         }
00405         
00406         // final step:  
00407         // and the tof and cer results together
00408         const UInt_t pid_result = (tof_pid&cer_pid);
00409         MSG("CalDetPID", Msg::kDebug)
00410              <<"result = "<<hex<<pid_result<<dec
00411              <<",  Identified as "
00412              <<CalDetParticleType::AsString(pid_result)<<endl;
00413         cpidh.SetPIDType(pid_result);
00414 
00415         MSG("CalDetPID", Msg::kDebug)
00416              <<"overlap_pid = "<<hex<<overlap_pid<<dec
00417              <<",  Non overlapping for "
00418              <<CalDetParticleType::AsString(overlap_pid)<<endl;
00419         
00420         cpidh.SetNoOverlapBits(overlap_pid);
00421 
00422         MSG("CalDetPID", Msg::kDebug)
00423              <<"ctw_pid = "<<hex<<ctw_pid<<dec
00424              <<",  Cer in time window for "
00425              <<CalDetParticleType::AsString(ctw_pid)<<endl;
00426 
00427         cpidh.SetInCERTimeBits(ctw_pid);
00428 
00429         Bool_t nooverlap;
00430         if((pid_result&overlap_pid)==pid_result) nooverlap=kTRUE;
00431         else nooverlap=kFALSE;
00432         MSG("CalDetPID", Msg::kDebug)
00433              <<"NoOverlap = "<<( (nooverlap) ? "true": "false")<<endl;
00434         cpidh.SetNoOverlap(nooverlap);
00435 
00436         Bool_t incertime;
00437         if((pid_result&ctw_pid)==pid_result) incertime=kTRUE;
00438         else incertime=kFALSE;
00439         MSG("CalDetPID", Msg::kDebug)
00440              <<"InCERTime = "<<( (incertime) ? "true": "false")<<endl;
00441         cpidh.SetInCERTime(incertime);
00442 
00443         cpidh.SetOLChi2(chi2);
00444    }
00445    else if(ismc){
00446         cpidh.SetPIDType(CalDetParticleType::kUnknown);
00447         cpidh.SetNoOverlapBits(0);
00448         cpidh.SetInCERTimeBits(0);
00449         cpidh.SetNoOverlap(kTRUE);
00450         cpidh.SetInCERTime(kTRUE);
00451         cpidh.SetOLChi2(0.);
00452         // look for the primary in stdhep
00453         const SimSnarlRecord* simrec = dynamic_cast<const SimSnarlRecord*>
00454             (cx.GetMom()->GetFragment("SimSnarlRecord"));
00455         if(simrec){
00456             const TClonesArray* stdhep = dynamic_cast<const TClonesArray*> 
00457                 (simrec->FindComponent("TClonesArray", "StdHep"));
00458             if(stdhep){
00459                 // take the first in the list
00460                 const TParticle* part = dynamic_cast<const TParticle*>
00461                     (stdhep->At(0));
00462                 if(part){
00463                     Int_t ipdg = part->GetPdgCode();
00464                     if(abs(ipdg)==11) 
00465                         cpidh.SetPIDType(CalDetParticleType::kElectron);
00466                     else if(abs(ipdg)==13) 
00467                         cpidh.SetPIDType(CalDetParticleType::kMuon);
00468                     else if(abs(ipdg)==211) 
00469                         cpidh.SetPIDType(CalDetParticleType::kPion);
00470                     else if(abs(ipdg)==321) 
00471                         cpidh.SetPIDType(CalDetParticleType::kKaon);
00472                     else if(abs(ipdg)==2212) 
00473                         cpidh.SetPIDType(CalDetParticleType::kProton);
00474                     }
00475                 }
00476             }
00477        }
00478 
00479 
00480 }

void AlgCalDetPID::Trace ( const char *  c  )  const [virtual]

Reimplemented from AlgBase.

Definition at line 483 of file AlgCalDetPID.cxx.

00484 {
00485 }


Member Data Documentation

Bool_t AlgCalDetPID::fFirstEvent [private]

Definition at line 45 of file AlgCalDetPID.h.

Referenced by RunAlg().

Definition at line 47 of file AlgCalDetPID.h.

Referenced by IsAnOverlap().

Definition at line 46 of file AlgCalDetPID.h.

Referenced by IsAnOverlap().

Int_t AlgCalDetPID::fNAllowedOut [private]

Definition at line 48 of file AlgCalDetPID.h.

Referenced by IsAnOverlap().

Int_t AlgCalDetPID::fNRequiredIn [private]

Definition at line 49 of file AlgCalDetPID.h.

Referenced by IsAnOverlap().

TH1F* AlgCalDetPID::fThist[2] [private]

Definition at line 52 of file AlgCalDetPID.h.

Referenced by HitTimeAnal(), OverlapLikelihood(), and RunAlg().

TFile* AlgCalDetPID::fTimeFile [private]

Definition at line 51 of file AlgCalDetPID.h.

Referenced by RunAlg(), and ~AlgCalDetPID().

Int_t AlgCalDetPID::fTimeFileNum [private]

Definition at line 50 of file AlgCalDetPID.h.

Referenced by RunAlg().


The documentation for this class was generated from the following files:

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1