TimingVarsAna Class Reference

#include <TimingVarsAna.h>

Inheritance diagram for TimingVarsAna:
NueAnaBase

List of all members.

Public Member Functions

 TimingVarsAna (TimingVars &tv)
virtual ~TimingVarsAna ()
void Analyze (int evtn, RecRecordImp< RecCandHeader > *srobj)

Public Attributes

TH1F * TimeHstShw
TH1F * TimeHstTrk
TH1F * TotalHst

Private Attributes

Detector::Detector_t fDetectorType
TimingVarsfTimingVars

Detailed Description

Definition at line 12 of file TimingVarsAna.h.


Constructor & Destructor Documentation

TimingVarsAna::TimingVarsAna ( TimingVars tv  ) 

Definition at line 28 of file TimingVarsAna.cxx.

References TimeHstShw, TimeHstTrk, and TotalHst.

00028                                           :
00029    fTimingVars(tv)
00030 {
00031 
00032 //    TimeHst = new THStack("TimeHist","Digit Times(ns) Red - Track, Blue - Shower.");
00033     TimeHstTrk = new TH1F();
00034     TimeHstTrk->SetName("TimeHstTrk");
00035     TimeHstShw = new TH1F();
00036     TimeHstShw->SetName("TimeHstShw");
00037 
00038     TotalHst = new TH1F();
00039     TotalHst->SetName("TotalHst");
00040 
00041 //    plot = new TCanvas("plot","plot",200,10,900,600);
00042 
00043 }

TimingVarsAna::~TimingVarsAna (  )  [virtual]

Definition at line 46 of file TimingVarsAna.cxx.

References TimeHstShw, TimeHstTrk, and TotalHst.

00047 {
00048   if(TimeHstTrk) delete TimeHstTrk;
00049   if(TimeHstShw) delete TimeHstShw;
00050   if(TotalHst) delete TotalHst;
00051 }


Member Function Documentation

void TimingVarsAna::Analyze ( int  evtn,
RecRecordImp< RecCandHeader > *  srobj 
) [virtual]

how many bins are over 40 ph

Implements NueAnaBase.

Definition at line 53 of file TimingVarsAna.cxx.

References UgliStripHandle::ClearFiber(), NtpSRTrack::ds, MuELoss::e, fDetectorType, fTimingVars, VldContext::GetDetector(), SntpHelpers::GetEvent(), UgliStripHandle::GetHalfLength(), RecRecordImp< T >::GetHeader(), SntpHelpers::GetShower(), SntpHelpers::GetShowerIndex(), SntpHelpers::GetSlice(), RecPhysicsHeader::GetSnarl(), VHS::GetStrip(), UgliGeomHandle::GetStripHandle(), SntpHelpers::GetTrack(), SntpHelpers::GetTrackIndex(), RecHeader::GetVldContext(), UgliStripHandle::GlobalToLocal(), ReleaseType::IsCedar(), ReleaseType::IsDogwood(), UgliGeomHandle::IsValid(), ReleaseType::kBirch, ReleaseType::kCedar, ReleaseType::kDogwood, Msg::kError, Detector::kFar, Detector::kNear, PlaneView::kU, PlaneView::kV, Msg::kWarning, StripEnd::kWest, MSG, NtpSRPlane::n, NtpSRTrack::nstrip, NtpSRShower::nstrip, NtpSRSlice::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRTrack::plane, NtpSRStrip::planeview, TimingVars::ShwBinsOver40, TimingVars::ShwbinsPrior, TimingVars::ShwEventPulseHeight, TimingVars::ShwMaxBin, TimingVars::ShwMaxBinCont, TimingVars::ShwPercentofTotal, TimingVars::ShwSecondMaxBin, NtpSRShower::stp, NtpSRSlice::stp, NtpSRTrack::stp, NtpSRTrack::stpds, NtpSRShower::stpt0, NtpSRTrack::stpt0, NtpSRShower::stpt1, NtpSRTrack::stpt1, NtpSRTrack::stpx, NtpSRTrack::stpy, NtpSRTrack::stpz, NtpSRStrip::strip, NtpSRStrip::time0, NtpSRStrip::time1, TimeHstShw, TimeHstTrk, TimingVars::TotalBinsOver40, TimingVars::TotalbinsPrior, TimingVars::TotalEventPulseHeight, TotalHst, TimingVars::TotalMaxBin, TimingVars::TotalMaxBinCont, TimingVars::TotalPercentofTotal, TimingVars::TotalSecondMaxBin, NtpSRStrip::tpos, TimingVars::TrkBinsOver40, TimingVars::TrkbinsPrior, TimingVars::TrkEventPulseHeight, TimingVars::TrkMaxBin, TimingVars::TrkMaxBinCont, TimingVars::TrkPercentofTotal, TimingVars::TrkSecondMaxBin, PropagationVelocity::Velocity(), UgliStripHandle::WlsPigtail(), and NtpSRStrip::z.

Referenced by NueRecordAna::Analyze().

00053                                                                        {
00054   if(srobj==0){
00055     return;
00056   }
00057   if(((dynamic_cast<NtpStRecord *>(srobj))==0)&&
00058      ((dynamic_cast<NtpSRRecord *>(srobj))==0)){
00059     return;
00060   }
00061 
00062   //Reset(srobj->GetHeader().GetSnarl(),evtn);
00063   NtpStRecord *st = dynamic_cast<NtpStRecord *>(srobj);
00064   NtpSREvent *event = SntpHelpers::GetEvent(evtn,srobj);
00065 
00066  
00067   if(!event){
00068       MSG("ShwfitAna",Msg::kError)<<"Couldn't get event "<<evtn
00069                                    <<" from Snarl "<<srobj->GetHeader().GetSnarl()<<endl;
00070       return;
00071    }
00072 
00073    VldContext vc=st->GetHeader().GetVldContext();
00074 
00075    //VldTimeStamp vldts = vc.GetTimeStamp();
00076    UgliGeomHandle ugh(vc);
00077    //UgliGeomHandle ugh = gMint->GetUgliGeomHandle();
00078    if (! ugh.IsValid()) {
00079      MSG("NueDisplayModule",Msg::kWarning) << "Got invalid Ugli\n";
00080      //return;
00081    }
00082 
00083   fDetectorType = vc.GetDetector();
00084   ReleaseType::Release_t fRel;
00085   string relName = st->GetTitle();
00086   string reco = relName.substr(0,relName.find_first_of("("));
00087   if(reco == "CEDAR") fRel = ReleaseType::kCedar;
00088   if(reco == "DOGWOOD") fRel = ReleaseType::kDogwood;
00089   else fRel = ReleaseType::kBirch;    
00090   int ntrks = event->ntrack;
00091       //timing histogram
00092     //find range
00093     double tmin=0, tmax=0;
00094     bool first = true;
00095     NtpSRSlice *slice = SntpHelpers::GetSlice(event->slc,st);
00096 
00097     //shamelessly stolen from Mad
00098     if (slice){//slice
00099       float highest_plane = 0;
00100       float lowest_plane = 500;
00101       float highest_strip0 = 0;
00102       float lowest_strip0 = 192;
00103       float highest_strip1 = 0;
00104       float lowest_strip1 = 192;
00105 
00106       float highest_z = 0.;
00107       float lowest_z = 30.;
00108       float highest_t0 = -4.0;
00109       float lowest_t0 = 4.0;
00110       float highest_t1 = -4.0;
00111       float lowest_t1 = 4.0;
00112 
00113       for (int i = 0; i<slice->nstrip; i++){//loop over strips
00114         NtpSRStrip *strip = SntpHelpers::GetStrip(slice->stp[i],srobj);
00115         if(strip == 0) continue;
00116         double t = strip->time0;
00117         if (t>-100){
00118           if (first){
00119             tmin = tmax = t;
00120             first = false;
00121           }
00122           if (t < tmin) tmin = t;
00123           if (t > tmax) tmax = t;
00124         }
00125         t = strip->time1;
00126         if (t>-100){
00127           if (first){
00128             tmin = tmax = t;
00129             first = false;
00130           }
00131           if (t < tmin) tmin = t;
00132           if (t > tmax) tmax = t;
00133         }
00134         int tempo_pln = strip->plane;
00135         int tempo_stp = strip->strip;
00136         float tempo_tpos = strip->tpos;
00137         if(tempo_pln<lowest_plane) {
00138           lowest_plane=tempo_pln;
00139           lowest_z=strip->z;
00140         }
00141         if(tempo_pln>highest_plane) {
00142           highest_plane=tempo_pln;
00143           highest_z=strip->z;
00144         }
00145 
00146         if (strip->planeview==PlaneView::kU){
00147          // fSlcUZ->Fill(strip->z,strip->tpos,(strip->ph0.sigcor + strip->ph1.sigcor)/SIGMAPMEU);
00148           if(tempo_tpos<lowest_t0) {
00149             lowest_strip0=tempo_stp;
00150             lowest_t0=tempo_tpos;
00151           }
00152           if(tempo_tpos>highest_t0) {
00153             highest_strip0=tempo_stp;
00154             highest_t0=tempo_tpos;
00155           }
00156         }
00157         else if (strip->planeview==PlaneView::kV){
00158           //fSlcVZ->Fill(strip->z,strip->tpos,(strip->ph0.sigcor + strip->ph1.sigcor)/SIGMAPMEU);
00159           if(tempo_tpos<lowest_t1) {
00160             lowest_strip1=tempo_stp;
00161             lowest_t1=tempo_tpos;
00162           }
00163           if(tempo_tpos>highest_t1) {
00164             highest_strip1=tempo_stp;
00165             highest_t1=tempo_tpos;
00166           }
00167         }
00168       }
00169       if(lowest_plane-10>=0) {
00170         lowest_plane-=10;
00171         lowest_z-=10.*0.06;
00172       }
00173       else {
00174         lowest_plane=0;
00175         lowest_z=0.;
00176       }
00177 
00178       if(lowest_strip0-5>=0) {
00179         lowest_strip0-=5;
00180         lowest_t0-=5.*0.041;
00181       }
00182       else {
00183         lowest_strip0=0;
00184         lowest_t0=-4.0;
00185       }
00186 
00187       if(lowest_strip1-5>=0) {
00188         lowest_strip1-=5;
00189         lowest_t1-=5.*0.041;
00190       }
00191       else {
00192         lowest_strip1=0;
00193         lowest_t1=-4.0;
00194       }
00195 
00196       if(highest_plane+10<=485) {
00197         highest_plane+=10;
00198         highest_z+=10.*0.06;
00199       }
00200       else {
00201         highest_plane=485;
00202         highest_z=30.;
00203       }
00204 
00205       if(highest_strip0+5<=191) {
00206         highest_strip0+=5;
00207         highest_t0+=5.*0.041;
00208       }
00209       else {
00210         highest_strip0=191;
00211         highest_t0=4.0;
00212       }
00213 
00214       if(highest_strip1+5<=191) {
00215         highest_strip1+=5;
00216         highest_t1+=5.*0.041;
00217       }
00218       else {
00219         highest_strip1=191;
00220         highest_t1=4.0;
00221       }
00222     }
00223     // give some buffer at either end...
00224     tmin -= 50e-9;
00225     tmax += 20e-9;
00226 
00227     double eps = 1.0e-8;
00228     if (tmin == tmax) { tmin -= eps; tmax += eps; }
00229 
00230     // for far detector, suppress display of pre-trigger time interval
00231     if (fDetectorType == Detector::kFar){
00232       if ((tmax - tmin)*1e9>1500) tmin = tmax - 1500e-9;
00233     }
00234 
00235     TimeHstTrk->SetBins(50,0,100);
00236     TimeHstShw->SetBins(50,0,100);
00237     TotalHst->SetBins(50,0,100);
00238     //record track hits
00239     //int ntrks = event->ntrack;
00240     if (ntrks){//if(ntrks)
00241       int trkidx = -1;
00242       int trkplanes = -1;
00243       for (int i = 0; i<ntrks; i++){
00244         int index = SntpHelpers::GetTrackIndex(i,event);
00245         NtpSRTrack *track = SntpHelpers::GetTrack(index,st);
00246         if(track == 0) continue;
00247         if (track->plane.n>trkplanes){
00248           trkplanes = track->plane.n;
00249           trkidx = index;
00250         }
00251         for (int istp = 0; istp<track->nstrip; istp++){
00252           //ftrkshw->Fill(track->stpx[istp],track->stpy[istp],1);
00253           NtpSRStrip *strip = SntpHelpers::GetStrip(track->stp[istp],st);
00254           if(strip == 0) continue;
00255           TimeHstTrk->Fill((track->stpt0[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph0.pe);
00256           TimeHstTrk->Fill((track->stpt1[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph1.pe);
00257           TotalHst->Fill((track->stpt0[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph0.pe);
00258           TotalHst->Fill((track->stpt1[istp]-tmin-(track->ds-track->stpds[istp])/3e8)/1e-9,strip->ph1.pe);
00259           //record track strip information
00260           if(strip->planeview==PlaneView::kU){
00261             //fTrkUZ->Fill(strip->z,strip->tpos,100000);//paranoia
00262           }
00263           else if(strip->planeview==PlaneView::kV){
00264             //fTrkVZ->Fill(strip->z,strip->tpos,100000);
00265           }
00266         }
00267       }
00268       vector<double> spathLength;
00269       vector<double> st0;
00270       if (trkidx != -1 && fDetectorType == Detector::kNear){
00271         NtpSRTrack *track = SntpHelpers::GetTrack(trkidx,st);
00272         for (int istp = 0; istp<track->nstrip; istp++){
00273           NtpSRStrip *strip = SntpHelpers::GetStrip(track->stp[istp],st);
00274           if(strip == 0) continue;
00275           spathLength.push_back(track->ds-track->stpds[istp]);
00276           PlexStripEndId seid(Detector::kNear,strip->plane,strip->strip,StripEnd::kWest);
00277           UgliStripHandle stripHandle = ugh.GetStripHandle(seid);
00278           float halfLength = stripHandle.GetHalfLength();
00279           const TVector3 ghitxyz(track->stpx[istp],track->stpy[istp],track->stpz[istp]);
00280           TVector3 lhitxyz = stripHandle.GlobalToLocal(ghitxyz);
00281           float fiberDist = (halfLength - lhitxyz.x() + stripHandle.ClearFiber(StripEnd::kWest) + stripHandle.WlsPigtail(StripEnd::kWest));
00282           //using strip time, I don't understand track time TJ
00283           //st0.push_back(track->stpt1[istp]-fiberDist/PropagationVelocity::Velocity());
00284           st0.push_back(strip->time1-fiberDist/PropagationVelocity::Velocity());
00285           //cout<<strip->plane<<" "<<strip->strip<<" "<<track->ds-track->stpds[istp]<<Form(" %.9f %.9f %.9f",track->stpt1[istp],strip->time1,strip->time1-fiberDist/PropagationVelocity::Velocity())<<endl;
00286         }
00287         float trms1 = 0;
00288         float trms2 = 0;
00289 
00290         trms1/=st0.size();
00291         trms2/=st0.size();
00292         trms1=sqrt(trms1)*1e9;
00293         trms2=sqrt(trms2)*1e9;
00294       }
00295     }
00296 
00297     //record shower hits
00298     int nshws = event->nshower;
00299     if (nshws){
00300       for (int i = 0; i<nshws; i++){
00301         int index = SntpHelpers::GetShowerIndex(i,event);
00302         NtpSRShower *shower = SntpHelpers::GetShower(index,st);
00303         if(shower == 0) continue;
00304         for (int istp = 0; istp<shower->nstrip; istp++){
00305           NtpSRStrip *strip = SntpHelpers::GetStrip(shower->stp[istp],st);
00306           if(strip == 0) continue;
00307           if(ReleaseType::IsCedar(fRel)||ReleaseType::IsDogwood(fRel)){
00308             TimeHstShw->Fill((shower->stpt0[istp]-tmin)/1e-9,strip->ph0.pe);
00309             TimeHstShw->Fill((shower->stpt1[istp]-tmin)/1e-9,strip->ph1.pe);
00310             TotalHst->Fill((shower->stpt0[istp]-tmin)/1e-9,strip->ph0.pe);
00311             TotalHst->Fill((shower->stpt1[istp]-tmin)/1e-9,strip->ph1.pe);
00312             //cout<<"Bin: "<<(shower->stpt1[istp]-tmin)/1e-9<<" Weight: "<<strip->ph0.pe+strip->ph1.pe<<endl;
00313           }else {
00314             TimeHstShw->Fill((strip->time0-tmin)/1e-9,strip->ph0.pe);
00315             TimeHstShw->Fill((strip->time1-tmin)/1e-9,strip->ph1.pe);
00316             TotalHst->Fill((strip->time0-tmin)/1e-9,strip->ph0.pe);
00317             TotalHst->Fill((strip->time1-tmin)/1e-9,strip->ph1.pe);
00318           }
00319           //record shower strip information
00320           if(strip->planeview==PlaneView::kU){
00321           //  fShwUZ->Fill(strip->z,strip->tpos,100000);//paranoia
00322           }
00323           else if(strip->planeview==PlaneView::kV){
00324            // fShwVZ->Fill(strip->z,strip->tpos,100000);
00325           }
00326         }
00327       }
00328     }
00329 //Max bins and contents
00330     Int_t ShwMaxBin = TimeHstShw->GetMaximumBin();
00331     Float_t ShwMaxBinCont = TimeHstShw->GetBinContent(ShwMaxBin);
00332 
00333     Int_t TrkMaxBin = TimeHstTrk->GetMaximumBin();
00334     Float_t TrkMaxBinCont = TimeHstTrk->GetBinContent(TrkMaxBin);
00335 
00336     Int_t TotalMaxBin = TotalHst->GetMaximumBin();
00337     Float_t TotalMaxBinCont = TotalHst->GetBinContent(TotalMaxBin);
00338 
00339 //Total Pulse height
00340     Float_t ShwEventPulseHeight = TimeHstShw->Integral();
00341     Float_t TrkEventPulseHeight = TimeHstTrk->Integral();
00342     Float_t TotalEventPulseHeight = TotalHst->Integral();
00343 
00344 //Second max bin contents
00345     Float_t ShwSecondMaxBin;
00346     if(TimeHstShw->GetBinContent(ShwMaxBin+1)>TimeHstShw->GetBinContent(ShwMaxBin-1)){
00347       ShwSecondMaxBin = TimeHstShw->GetBinContent(ShwMaxBin+1);
00348     }else{
00349       ShwSecondMaxBin = TimeHstShw->GetBinContent(ShwMaxBin-1);
00350     }
00351     Float_t TrkSecondMaxBin;
00352     if(TimeHstTrk->GetBinContent(TrkMaxBin+1)>TimeHstTrk->GetBinContent(TrkMaxBin-1)){
00353       TrkSecondMaxBin = TimeHstTrk->GetBinContent(TrkMaxBin+1);
00354     }else{
00355       TrkSecondMaxBin = TimeHstTrk->GetBinContent(TrkMaxBin-1);
00356     }
00357     Float_t TotalSecondMaxBin;
00358     if(TotalHst->GetBinContent(TotalMaxBin+1)>TotalHst->GetBinContent(TotalMaxBin-1)){
00359       TotalSecondMaxBin = TotalHst->GetBinContent(TotalMaxBin+1);
00360     }else{
00361       TotalSecondMaxBin = TotalHst->GetBinContent(TotalMaxBin-1);
00362     }
00363 
00364 //Percent of pulse height in Max 2 bins
00365     Float_t TrkPercentofTotal = -9999.99;
00366     Float_t ShwPercentofTotal = -9999.99;
00367     Float_t TotalPercentofTotal = -9999.99;
00368 
00369     if (ShwEventPulseHeight>0){
00370       ShwPercentofTotal = (ShwSecondMaxBin + ShwMaxBinCont)/(ShwEventPulseHeight);
00371     }
00372     if (TrkEventPulseHeight>0){
00373       TrkPercentofTotal = (TrkSecondMaxBin + TrkMaxBinCont)/(TrkEventPulseHeight);
00374     }
00375     if (TotalEventPulseHeight>0){
00376       TotalPercentofTotal = (TotalSecondMaxBin + TotalMaxBinCont)/(TotalEventPulseHeight);
00377     }
00378 
00379 //how many of the two bins prior to the max bin are at least 20% of the max bin ph
00380     Int_t ShwbinsPrior = 0;
00381     if(TimeHstShw->GetBinContent(ShwMaxBin-1)>(0.2*ShwMaxBinCont)){ShwbinsPrior++;}
00382     if(TimeHstShw->GetBinContent(ShwMaxBin-2)>(0.2*ShwMaxBinCont)){ShwbinsPrior++;}
00383     Int_t TrkbinsPrior = 0;
00384     if(TimeHstTrk->GetBinContent(TrkMaxBin-1)>(0.2*TrkMaxBinCont)){TrkbinsPrior++;}
00385     if(TimeHstTrk->GetBinContent(TrkMaxBin-2)>(0.2*TrkMaxBinCont)){TrkbinsPrior++;}
00386     Int_t TotalbinsPrior = 0;
00387     if(TotalHst->GetBinContent(TotalMaxBin-1)>(0.2*TotalMaxBinCont)){TotalbinsPrior++;}
00388     if(TotalHst->GetBinContent(TotalMaxBin-2)>(0.2*TotalMaxBinCont)){TotalbinsPrior++;}
00389     
00391     Int_t ShwBinsOver40 = 0;
00392     Int_t TrkBinsOver40 = 0;
00393     Int_t TotalBinsOver40 = 0;
00394     for (int k=0;k<50;k++){
00395       if(TimeHstShw->GetBinContent(k)>20)ShwBinsOver40++;
00396       if(TimeHstTrk->GetBinContent(k)>20)TrkBinsOver40++;
00397       if(TotalHst->GetBinContent(k)>20)TotalBinsOver40++;
00398     }
00399 
00400     fTimingVars.ShwMaxBin = ShwMaxBin;
00401     fTimingVars.TrkMaxBin = TrkMaxBin;
00402     fTimingVars.TotalMaxBin = TotalMaxBin;
00403 
00404     fTimingVars.ShwMaxBinCont = ShwMaxBinCont;
00405     fTimingVars.TrkMaxBinCont = TrkMaxBinCont;
00406     fTimingVars.TotalMaxBinCont = TotalMaxBinCont;
00407 
00408     fTimingVars.ShwSecondMaxBin = ShwSecondMaxBin;
00409     fTimingVars.TrkSecondMaxBin = TrkSecondMaxBin;
00410     fTimingVars.TotalSecondMaxBin = TotalSecondMaxBin;
00411 
00412     fTimingVars.ShwEventPulseHeight = ShwEventPulseHeight;
00413     fTimingVars.TrkEventPulseHeight = TrkEventPulseHeight;
00414     fTimingVars.TotalEventPulseHeight = TotalEventPulseHeight;
00415 
00416     fTimingVars.ShwPercentofTotal = ShwPercentofTotal;
00417     fTimingVars.TrkPercentofTotal = TrkPercentofTotal;
00418     fTimingVars.TotalPercentofTotal = TotalPercentofTotal;
00419 
00420     fTimingVars.ShwbinsPrior = ShwbinsPrior;
00421     fTimingVars.TrkbinsPrior = TrkbinsPrior;
00422     fTimingVars.TotalbinsPrior = TotalbinsPrior;
00423 
00424     fTimingVars.ShwBinsOver40 = ShwBinsOver40;
00425     fTimingVars.TrkBinsOver40 = TrkBinsOver40;
00426     fTimingVars.TotalBinsOver40 = TotalBinsOver40;
00427     /*
00428     char name[100];
00429     sprintf(name,"/home/danche/ViewScan5/plot/ShwHst_snarl_%d.png",st->GetHeader().GetSnarl());
00430 
00431     TimeHstShw->Draw();
00432     plot->Print(name);
00433     */
00434 }


Member Data Documentation

Definition at line 28 of file TimingVarsAna.h.

Referenced by Analyze().

Definition at line 29 of file TimingVarsAna.h.

Referenced by Analyze().

Definition at line 21 of file TimingVarsAna.h.

Referenced by Analyze(), TimingVarsAna(), and ~TimingVarsAna().

Definition at line 22 of file TimingVarsAna.h.

Referenced by Analyze(), TimingVarsAna(), and ~TimingVarsAna().

Definition at line 23 of file TimingVarsAna.h.

Referenced by Analyze(), TimingVarsAna(), and ~TimingVarsAna().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1