#include <TimingVarsAna.h>
Inheritance diagram for TimingVarsAna:

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 |
| TimingVars & | fTimingVars |
|
|
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 }
|
|
|
Definition at line 46 of file TimingVarsAna.cxx. 00047 {
00048 if(TimeHstTrk) delete TimeHstTrk;
00049 if(TimeHstShw) delete TimeHstShw;
00050 if(TotalHst) delete TotalHst;
00051 }
|
|
||||||||||||
|
how many bins are over 40 ph Implements NueAnaBase. Definition at line 53 of file TimingVarsAna.cxx. References UgliStripHandle::ClearFiber(), NtpSRTrack::ds, fDetectorType, fTimingVars, VldContext::GetDetector(), SntpHelpers::GetEvent(), UgliStripHandle::GetHalfLength(), RecRecordImp< T >::GetHeader(), SntpHelpers::GetShower(), SntpHelpers::GetShowerIndex(), SntpHelpers::GetSlice(), RecPhysicsHeader::GetSnarl(), SntpHelpers::GetStrip(), UgliGeomHandle::GetStripHandle(), SntpHelpers::GetTrack(), SntpHelpers::GetTrackIndex(), RecHeader::GetVldContext(), UgliStripHandle::GlobalToLocal(), ReleaseType::IsCedar(), ReleaseType::IsDogwood(), UgliGeomHandle::IsValid(), MSG, NtpSRPlane::n, NtpSREvent::nshower, NtpSRShower::nstrip, NtpSRTrack::nstrip, NtpSRSlice::nstrip, NtpSREvent::ntrack, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRTrack::plane, NtpSRStrip::plane, NtpSRStrip::planeview, TimingVars::ShwBinsOver40, TimingVars::ShwbinsPrior, TimingVars::ShwEventPulseHeight, TimingVars::ShwMaxBin, TimingVars::ShwMaxBinCont, TimingVars::ShwPercentofTotal, TimingVars::ShwSecondMaxBin, NtpSREvent::slc, NtpSRShower::stp, NtpSRTrack::stp, NtpSRSlice::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, 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 }
|
|
|
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(), and TimingVarsAna(). |
|
|
Definition at line 22 of file TimingVarsAna.h. Referenced by Analyze(), and TimingVarsAna(). |
|
|
Definition at line 23 of file TimingVarsAna.h. Referenced by Analyze(), and TimingVarsAna(). |
1.3.9.1